1 | //last change: 13.02.2001 (Eric Westenberger) |
---|
2 | /////////////////////////////////////////////////////////////////////////////// |
---|
3 | version="$Id: graphics.lib,v 1.13 2009-04-06 09:17:01 seelisch Exp $"; |
---|
4 | category="Visualization"; |
---|
5 | info=" |
---|
6 | LIBRARY: graphics.lib Procedures for graphical output using Mathematica |
---|
7 | AUTHOR: Christian Gorzel, gorzelc@math.uni-muenster.de |
---|
8 | |
---|
9 | PROCEDURES: |
---|
10 | staircase(fname,I); Mathematica text for displaying staircase of I |
---|
11 | mathinit(); string for loading Mathematica's ImplicitPlot |
---|
12 | mplot(fname,I[# s]); Mathematica text for various plots |
---|
13 | "; |
---|
14 | |
---|
15 | /////////////////////////////////////////////////////////////////////////////// |
---|
16 | |
---|
17 | proc staircase(ideal I) |
---|
18 | "USAGE: staircase(I); I an ideal in two variables |
---|
19 | RETURN: string with Mathematica input for displaying staircase diagrams of an |
---|
20 | ideal I, i.e. exponent vectors of the initial ideal of I |
---|
21 | NOTE: ideal I should be given by a standard basis. Copy and |
---|
22 | paste the result into a Mathematica notebook. |
---|
23 | EXAMPLE: example staircase; shows an example |
---|
24 | " |
---|
25 | { |
---|
26 | intvec v; |
---|
27 | int maxx, maxy; |
---|
28 | list l; |
---|
29 | string s; |
---|
30 | string el; |
---|
31 | |
---|
32 | if(nvars(basering)!=2) |
---|
33 | { "-- Error: need two variables "; |
---|
34 | return(); |
---|
35 | } |
---|
36 | if (not(attrib(I,"isSB"))) |
---|
37 | { " -- Warning: Ideal should be a standardbasis "; newline; } |
---|
38 | |
---|
39 | for(int i=1; i<=ncols(I); i++) |
---|
40 | { |
---|
41 | if (i>1) { el = el + ",";} |
---|
42 | v = leadexp(I[i]); |
---|
43 | if (v[1] > maxx) { maxx = v[1];} |
---|
44 | if (v[2] > maxy) { maxy = v[2];} |
---|
45 | el = el + "{" + string(v) + "}"; |
---|
46 | } |
---|
47 | el = "{" + el + "}"; |
---|
48 | maxx = maxx + 3; |
---|
49 | maxy = maxy + 3; |
---|
50 | |
---|
51 | s = newline + |
---|
52 | "Show[Graphics[{" + newline + |
---|
53 | "{GrayLevel[0.5],Map[Rectangle[#,{" + |
---|
54 | string(maxx) + "," + string(maxy) + "}] &, " + el + "]}," + newline + |
---|
55 | "{PointSize[0.03], Map[Point," + el + "]}," + newline + |
---|
56 | "Table[Circle[{i,j},0.1],{i,0," + |
---|
57 | string(maxx) + "},{j,0," + string(maxy) + "}]}," + newline + |
---|
58 | " Axes->True,AspectRatio->Automatic]]"; |
---|
59 | s = s + endstr(fname); |
---|
60 | return(s); |
---|
61 | } |
---|
62 | example |
---|
63 | { "EXAMPLE:"; echo =2; |
---|
64 | ring r0 = 0,(x,y),ls; |
---|
65 | ideal I = -1x2y6-1x4y2, 7x6y5+1/2x7y4+6x4y6; |
---|
66 | staircase(std(I)); |
---|
67 | |
---|
68 | ring r1 = 0,(x,y),dp; |
---|
69 | ideal I = fetch(r0,I); |
---|
70 | staircase(std(I)); |
---|
71 | |
---|
72 | ring r2 = 0,(x,y),wp(2,3); |
---|
73 | ideal I = fetch(r0,I); |
---|
74 | staircase(std(I)); |
---|
75 | |
---|
76 | // Paste the output into a Mathematica notebook |
---|
77 | // active evalutation of the cell with SHIFT RETURN |
---|
78 | } |
---|
79 | /////////////////////////////////////////////////////////////////////////////// |
---|
80 | |
---|
81 | proc mathinit() |
---|
82 | "USAGE: mathinit(); |
---|
83 | RETURN: initializing string for loading Mathematica's ImplicitPlot |
---|
84 | EXAMPLE: example mathinit; shows an example |
---|
85 | " |
---|
86 | { |
---|
87 | // write("init.m","<< Graphics`ImplicitPlot`"); |
---|
88 | return("<< Graphics`ImplicitPlot`"); |
---|
89 | } |
---|
90 | example |
---|
91 | { "EXAMPLE:"; echo =2; |
---|
92 | mathinit(); |
---|
93 | |
---|
94 | // Paste the output into a Mathematica notebook |
---|
95 | // active evalutation of the cell with SHIFT RETURN |
---|
96 | } |
---|
97 | /////////////////////////////////////////////////////////////////////////////// |
---|
98 | |
---|
99 | static proc checkshort() |
---|
100 | { |
---|
101 | ring @r; |
---|
102 | } |
---|
103 | static proc determvars(ideal I) |
---|
104 | // determine the variables which are in the ideal I |
---|
105 | { |
---|
106 | intvec v; |
---|
107 | int i,j,k; |
---|
108 | |
---|
109 | k=1; |
---|
110 | for(j=1;j<=size(I);j++) |
---|
111 | { for(i=1;i<=nvars(basering);i++) |
---|
112 | { if(I[j]!=subst(I[j],var(i),0)) {v[k] = i; k++;} |
---|
113 | } |
---|
114 | } |
---|
115 | ring @r=0,x,ls; |
---|
116 | poly f; |
---|
117 | for(i=1;i<=size(v);i++) // sort VARS |
---|
118 | { f = f + x^v[i]; } |
---|
119 | v=0; |
---|
120 | for (i=1;i<=size(f);i++) |
---|
121 | {v[i]=leadexp(f[i])[1];} |
---|
122 | return(v); |
---|
123 | } |
---|
124 | /////////////////////////////////////////////////////////////////////////////// |
---|
125 | |
---|
126 | static proc endstr(string filename) |
---|
127 | { int i; |
---|
128 | string suffix,cmd,name; |
---|
129 | |
---|
130 | if(size(filename)) |
---|
131 | { |
---|
132 | for (i=size(filename);i;i--) |
---|
133 | { if (filename[i] == ".") {break;} |
---|
134 | } |
---|
135 | |
---|
136 | if (i>0) |
---|
137 | { suffix = filename[i,size(filename)-i+1]; |
---|
138 | name = ">" + filename[1,i-1]+ ".m"; |
---|
139 | } |
---|
140 | else { print("--Error: Suffix of filename incorrect"); return("");} |
---|
141 | // if (suffix ==".m") { cmd = "Display[\" " + filename + "\",% ]";} |
---|
142 | if (suffix ==".mps") { cmd = "Display[\" " + filename + "\",%] ";} |
---|
143 | if (suffix ==".ps") { cmd = "Display[\" ! psfix > " + filename + "\", %]";} |
---|
144 | if (suffix ==".eps") |
---|
145 | { cmd = "Display[\" ! psfix -epsf > " + filename + "\", %]";} |
---|
146 | |
---|
147 | } |
---|
148 | return(newline + cmd); |
---|
149 | } |
---|
150 | |
---|
151 | /////////////////////////////////////////////////////////////////////////////// |
---|
152 | proc mplot(string fname,ideal I,list #) |
---|
153 | "USAGE: mplot(fname, I [,I1,I2,..,s] ); fname=string; I,I1,I2,..=ideals, |
---|
154 | s=string representing the plot region.@* |
---|
155 | Use the ideals I1,I2,.. in order to produce multiple plots (they need |
---|
156 | to have the same number of entries as I!). |
---|
157 | RETURN: string, text with Mathematica commands to display a plot |
---|
158 | NOTE: The plotregion is defaulted to -1,1 around zero. |
---|
159 | For implicit given curves enter first the string returned by |
---|
160 | procedure mathinit into Mathematica in order to load ImplicitPlot. |
---|
161 | The following conventions for I are used: |
---|
162 | @format |
---|
163 | - ideal with 2 entries in one variable means a parametrised plane curve, |
---|
164 | - ideal with 3 entries in one variable means a parametrised space curve, |
---|
165 | - ideal with 3 entries in two variables means a parametrised surface, |
---|
166 | - ideal with 2 entries in two variables means an implicit curve |
---|
167 | given as I[1]==I[2], |
---|
168 | - ideal with 1 entry (or one polynomial) in two variables means |
---|
169 | an implicit curve given as f == 0, |
---|
170 | @end format |
---|
171 | EXAMPLE: example mplot; shows an example |
---|
172 | " |
---|
173 | { |
---|
174 | int i,j,k,mapping; |
---|
175 | int planecurve,spacecurve,implcrv,surface; |
---|
176 | intvec VARS,v; |
---|
177 | string xpart,ypart,zpart = "-1,1","-1,1","All"; |
---|
178 | string pstring,actstring,xname,yname,str,closing; |
---|
179 | string basr = nameof(basering); |
---|
180 | ideal J; |
---|
181 | |
---|
182 | if (ncols(I)>3) |
---|
183 | { "-- Error: can only draw upto dimension 3"; |
---|
184 | return(""); |
---|
185 | } |
---|
186 | ring @r = 0,(s,t),lp; |
---|
187 | ideal @J,@I; |
---|
188 | |
---|
189 | setring(`basr`); |
---|
190 | // def d = basering; |
---|
191 | // d; |
---|
192 | // listvar(d); |
---|
193 | |
---|
194 | str = newline; |
---|
195 | |
---|
196 | VARS = determvars(I); |
---|
197 | // "VARS: ";VARS; |
---|
198 | |
---|
199 | if (size(VARS)>2 or VARS==0) |
---|
200 | { "-- Error: Need some variables, but can only draw in 2 or 3 dimensions"; |
---|
201 | return(""); |
---|
202 | } |
---|
203 | if (size(matrix(I))==1 and size(VARS)==2) |
---|
204 | { i =size(I[1]); |
---|
205 | //I[2]=I[1][(i/ 2 + 1)..i]; I[2]; |
---|
206 | // I[1]=I[1][1..(i/ 2)]; I[1]; |
---|
207 | I[2]=0; |
---|
208 | } |
---|
209 | if (size(matrix(I))==2) |
---|
210 | { if (size(VARS)==1) {planecurve=1; str = str + "ParametricPlot[";} |
---|
211 | if (size(VARS)==2) {implcrv=1; str = str + "ImplicitPlot[";} |
---|
212 | } |
---|
213 | if (size(matrix(I))==3) |
---|
214 | { if (size(VARS)==1) {spacecurve=1;} |
---|
215 | if (size(VARS)==2) {surface=1;} |
---|
216 | str = str + "ParametricPlot3D["; |
---|
217 | } |
---|
218 | |
---|
219 | short = 0; |
---|
220 | |
---|
221 | pstring = string(I); |
---|
222 | |
---|
223 | // switch to another ring if necessary |
---|
224 | |
---|
225 | checkshort(); |
---|
226 | // "short: "; short; |
---|
227 | |
---|
228 | if (short!=1) // construct a map |
---|
229 | { |
---|
230 | mapping = 1; |
---|
231 | setring @r; |
---|
232 | @J = 0; |
---|
233 | for(i=1;i<=size(VARS);i++) |
---|
234 | { @J[VARS[i]]=var(i);} |
---|
235 | map phi = (`basr`,@J); |
---|
236 | @I = phi(I); |
---|
237 | short =0; |
---|
238 | pstring = string(@I); |
---|
239 | setring `basr`; |
---|
240 | } |
---|
241 | |
---|
242 | i = find(pstring,newline); |
---|
243 | while(i) |
---|
244 | {pstring[i]=" "; |
---|
245 | i = find(pstring,newline,i); |
---|
246 | } |
---|
247 | if (implcrv) |
---|
248 | { i = find(pstring,","); |
---|
249 | pstring = pstring[1,i-1] + "==" + pstring[i+1,size(pstring)-i]; |
---|
250 | } |
---|
251 | else |
---|
252 | { pstring = "{" + pstring + "}";} |
---|
253 | // "mapping "; mapping; |
---|
254 | if (mapping) |
---|
255 | { xname = "s"; |
---|
256 | if (size(VARS)==2) {yname="t";} |
---|
257 | } |
---|
258 | else |
---|
259 | { xname = varstr(VARS[1]); |
---|
260 | if (size(VARS)==2) {yname=varstr(VARS[2]);} |
---|
261 | } |
---|
262 | |
---|
263 | j =1; |
---|
264 | |
---|
265 | for(k=1;k<=size(#);k++) |
---|
266 | { if (typeof(#[k])=="ideal" or typeof(#[k])=="poly") |
---|
267 | { //#[k] = ideal(#[k]); |
---|
268 | v = determvars(#[k]); |
---|
269 | J = #[k]; |
---|
270 | short =0; |
---|
271 | if (size(matrix(J))==1 and size(VARS)==2 and implcrv) |
---|
272 | { i =size(J[1]); |
---|
273 | // J[2]=J[1][(i/ 2 + 1)..i]; |
---|
274 | // J[1]=J[1][1..(i/ 2)]; |
---|
275 | J[2] =0; |
---|
276 | } |
---|
277 | if ((v!= VARS) or (size(J)!=size(I))) |
---|
278 | { print("--Error: ---- "); |
---|
279 | return(); |
---|
280 | } |
---|
281 | else |
---|
282 | { if (mapping) |
---|
283 | { setring @r; |
---|
284 | short =0; |
---|
285 | actstring = string(phi(J)); |
---|
286 | setring(`basr`); |
---|
287 | } |
---|
288 | else {actstring = string(J);} |
---|
289 | i = find(actstring,newline); |
---|
290 | while(i) |
---|
291 | { actstring[i]=" "; |
---|
292 | i = find(actstring,newline,i); |
---|
293 | } |
---|
294 | if (implcrv) |
---|
295 | {i = find(actstring,","); |
---|
296 | actstring = actstring[1,i-1]+ "==" + actstring[i+1,size(actstring)-i]; |
---|
297 | pstring = pstring + "," + actstring; |
---|
298 | } |
---|
299 | else |
---|
300 | { pstring = pstring + ",{" + actstring +"}"; |
---|
301 | } |
---|
302 | |
---|
303 | } |
---|
304 | } |
---|
305 | if (typeof(#[k])=="string") |
---|
306 | { if (j==3) {zpart = #[k];j++;} |
---|
307 | if (j==2) {ypart = #[k];j++;} |
---|
308 | if (j==1) {xpart = #[k];j++;} |
---|
309 | } |
---|
310 | } |
---|
311 | |
---|
312 | if (spacecurve or planecurve or implcrv) |
---|
313 | { str = str + "{" + pstring + "},{" + xname + "," + xpart + "}";} |
---|
314 | if (implcrv and (j==3)) {str = str + ",{" + yname + "," + ypart + "}";} |
---|
315 | if (surface) |
---|
316 | { str = str + "{" + pstring + "},{" + xname + "," + xpart + "},{" |
---|
317 | + yname + "," + ypart + "}";} |
---|
318 | |
---|
319 | if (planecurve) {closing = "," + newline +" AspectRatio->Automatic";} |
---|
320 | if (implcrv) {closing = "," + newline + |
---|
321 | " AxesLabel->{\"" + varstr(VARS[1]) + "\",\"" + varstr(VARS[2]) + "\"}";} |
---|
322 | if (spacecurve) { closing = "," + newline + " ViewPoint->{1.3,-2.4,2}";} |
---|
323 | if (surface) |
---|
324 | {closing = "," +newline + |
---|
325 | " Boxed->True, Axes->True, ViewPoint->{1.3,-2.4,2}";} |
---|
326 | |
---|
327 | str = str + closing + "];" + endstr(fname); |
---|
328 | |
---|
329 | return(str); |
---|
330 | } |
---|
331 | example |
---|
332 | { "EXAMPLE:"; echo =2; |
---|
333 | // --------- plane curves ------------ |
---|
334 | ring rr0 = 0,x,dp; export rr0; |
---|
335 | |
---|
336 | ideal I = x3 + x, x2; |
---|
337 | ideal J = x2, -x+x3; |
---|
338 | mplot("",I,J,"-2,2"); |
---|
339 | |
---|
340 | // Paste the output into a Mathematica notebook |
---|
341 | // active evalutation of the cell with SHIFT RETURN |
---|
342 | |
---|
343 | pause("Hit RETURN to continue"); |
---|
344 | // --------- space curves -------------- |
---|
345 | I = x3,-1/10x3+x2,x2; |
---|
346 | mplot("",I); |
---|
347 | |
---|
348 | // Paste the output into a Mathematica notebook |
---|
349 | // active evalutation of the cell with SHIFT RETURN |
---|
350 | |
---|
351 | pause("Hit RETURN to continue"); |
---|
352 | // ----------- surfaces ------------------- |
---|
353 | ring rr1 = 0,(x,y),dp; export rr1; |
---|
354 | ideal J = xy,y,x2; |
---|
355 | mplot("",J,"-2,1","1,2"); |
---|
356 | |
---|
357 | // Paste the output into a Mathematica notebook |
---|
358 | // active evalutation of the cell with SHIFT RETURN |
---|
359 | kill rr0,rr1; |
---|
360 | } |
---|
361 | /////////////////////////////////////////////////////////////////////////////// |
---|