source: git/Singular/LIB/graphics.lib @ c7c9cb

spielwiese
Last change on this file since c7c9cb was c7c9cb, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: initial addition of "graphics.lib" git-svn-id: file:///usr/local/Singular/svn/trunk@1754 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.6 KB
Line 
1// $Id: graphics.lib,v 1.1 1998-05-13 17:26:11 Singular Exp $
2//
3// author : Christian Gorzel email: gorzelc@math.uni-muenster.de
4//
5///////////////////////////////////////////////////////////////////////////////
6version="$Id: graphics.lib,v 1.1 1998-05-13 17:26:11 Singular Exp $";
7info="
8LIBRARY: graphics.lib    PROCEDURES FOR GRAPHICS WITH MATHEMATICA
9
10 staircase(ideal I);  string of Mathematica input for displaying staircase of I
11 mathinit();          string for loading Mathematica's ImplicitPlot
12 plot(ideal I,list #);string with Mathematica text for various plots
13";
14
15///////////////////////////////////////////////////////////////////////////////
16
17proc staircase(ideal I)
18"USAGE:   staircase(); I ideal in two variables
19RETURN:  string with Mathematica input for displaying staircase (exponent
20         vectors of initial ideal) of I
21NOTE:    ideal I should be given by a standard basis
22         Copy and paste works only in a Mathematica notebook
23EXAMPLE: 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
60 return(s);
61}
62example
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
81proc mathinit()
82"USAGE:   mathinit();
83RETURN:  initializing string for loading Mathematica's ImplicitPlot
84EXAMPLE: example mathinit; shows an example
85"
86{
87  // write("init.m","<< Graphics`ImplicitPlot`");
88  return("<< Graphics`ImplicitPlot`");
89}
90example
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///////////////////////////////////////////////////////////////////////////////
98static proc checkshort()
99{
100  ring @r;
101}
102static proc determvars(ideal I)
103// determine the variables which are in the ideal I
104{
105  intvec v;
106  int i,j,k;
107
108  k=1;
109  for(j=1;j<=size(I);j++)
110  { for(i=1;i<=nvars(basering);i++)
111    { if(I[j]!=subst(I[j],var(i),0)) {v[k] = i; k++;}
112    }
113  }
114  ring @r=0,x,ls;
115  poly f;
116  for(i=1;i<=size(v);i++)     // sort VARS
117  { f = f + x^v[i]; }
118  v=0;
119  for (i=1;i<=size(f);i++)
120  {v[i]=leadexp(f[i])[1];}
121 return(v);
122}
123///////////////////////////////////////////////////////////////////////////////
124proc  plot(ideal I,list #)
125"USAGE:   plot(I[,#]); I ideal , #  strings representing the plot region or
126                                   further ideals
127RETURN:  string, text with Mathematica commands to display a plot
128NOTE:    The plotregion is defaulted to -1,1 around zero
129         (mostly interested in local situation in singularity theory)
130         For implicit given curves enter first the string returned by
131         proc mathinit into Mathematica in order to load ImplicitPlot
132
133         ideal with 2 entries in one variable means a parametrised plane curve
134         ideal with 3 entries in one variable means a parametrised space curve
135         ideal with 3 entries in two variables means a parametrised surface
136
137         ideal I with 2 generators in two variables means an implicit curve
138          given as I[1]==I[2]
139         ideal with 1 generator (or one polynomial) in two variables means
140          an implicit curve given as  f == 0
141
142EXAMPLE: example plot; shows an example
143"
144{
145  int i,j,k,mapping;
146  int planecurve,spacecurve,implcrv,surface;
147  intvec VARS,v;
148  string xpart,ypart,zpart = "-1,1","-1,1","All";
149  string pstring,actstring,xname,yname,str,closing;
150  string basr = nameof(basering);
151  ideal J;
152
153  ring @r = 0,(s,t),lp;
154  ideal @J,@I;
155
156  setring(`basr`);
157  // def d = basering;
158  // d;
159  // listvar(d);
160
161  str = newline;
162
163  VARS = determvars(I);
164   // "VARS: ";VARS;
165
166  if (size(matrix(I))==1 and size(VARS)==2)
167  { i =size(I[1]);
168   //I[2]=I[1][(i/ 2 + 1)..i]; I[2];
169   // I[1]=I[1][1..(i/ 2)]; I[1];
170   I[2]=0;
171  }
172  if (size(matrix(I))==2)
173  { if (size(VARS)==1) {planecurve=1; str = str + "ParametricPlot[";}
174    if (size(VARS)==2) {implcrv=1; str = str + "ImplicitPlot[";}
175  }
176  if (size(matrix(I))==3)
177  { if (size(VARS)==1) {spacecurve=1;}
178    if (size(VARS)==2) {surface=1;}
179    str =  str + "ParametricPlot3D[";
180  }
181
182  short = 0;
183
184  pstring = string(I);
185
186 // switch to another ring if necessary
187
188  checkshort();
189//  "short: "; short;
190
191  if (short!=1)      // construct a map
192  {
193    mapping = 1;
194    setring @r;
195    @J = 0;
196    for(i=1;i<=size(VARS);i++)
197    { @J[VARS[i]]=var(i);}
198    map phi = (`basr`,@J);
199    @I = phi(I);
200    short =0;
201    pstring = string(@I);
202    setring `basr`;
203  }
204
205  i = find(pstring,newline);
206  while(i)
207  {pstring[i]=" ";
208   i = find(pstring,newline,i);
209  }
210  if (implcrv)
211  { i = find(pstring,",");
212    pstring = pstring[1,i-1] + "==" + pstring[i+1,size(pstring)-i];
213  }
214  else
215  { pstring = "{" + pstring + "}";}
216//  "mapping "; mapping;
217  if (mapping)
218  { xname = "s";
219    if (size(VARS)==2) {yname="t";}
220  }
221  else
222  { xname = varstr(VARS[1]);
223    if (size(VARS)==2) {yname=varstr(VARS[2]);}
224  }
225
226  j =1;
227
228  for(k=1;k<=size(#);k++)
229  { if (typeof(#[k])=="ideal" or typeof(#[k])=="poly")
230    { //#[k] = ideal(#[k]);
231      v = determvars(#[k]);
232      J = #[k];
233      short =0;
234      if (size(matrix(J))==1 and size(VARS)==2 and implcrv)
235      { i =size(J[1]);
236      //  J[2]=J[1][(i/ 2 + 1)..i];
237      //  J[1]=J[1][1..(i/ 2)];
238        J[2] =0;
239      }
240      if ((v!= VARS) or (size(J)!=size(I)))
241      { print("--Error: ---- ");
242        return();
243      }
244      else
245      { if (mapping)
246        { setring @r;
247          short =0;
248          actstring = string(phi(J));
249          setring(`basr`);
250        }
251        else {actstring = string(J);}
252        i = find(actstring,newline);
253        while(i)
254        { actstring[i]=" ";
255         i = find(actstring,newline,i);
256        }
257        if (implcrv)
258        {i = find(actstring,",");
259         actstring = actstring[1,i-1]+ "==" + actstring[i+1,size(actstring)-i];
260         pstring = pstring + "," + actstring;
261        }
262        else
263        { pstring = pstring + ",{" + actstring +"}";
264        }
265
266      }
267    }
268    if (typeof(#[k])=="string")
269    {  if (j==3) {zpart = #[k];j++;}
270       if (j==2) {ypart = #[k];j++;}
271       if (j==1) {xpart = #[k];j++;}
272   }
273  }
274
275 if (spacecurve or planecurve or implcrv)
276 { str = str + "{" + pstring + "},{" + xname + "," + xpart + "}";}
277 if (implcrv and (j==3)) {str = str + ",{" + yname + "," + ypart + "}";}
278 if (surface)
279 { str = str + "{" + pstring + "},{" + xname + "," + xpart + "},{"
280                                     + yname + "," + ypart + "}";}
281
282  if (planecurve) {closing = "," + newline +" AspectRatio->Automatic";}
283  if (implcrv) {closing = "," + newline +
284   " AxesLabel->{\"" + varstr(VARS[1]) + "\",\"" + varstr(VARS[2]) + "\"}";}
285  if (spacecurve) { closing = "," + newline + " ViewPoint->{1.3,-2.4,2}";}
286  if (surface)
287  {closing = "," +newline +
288              " Boxed->True, Axes->True, ViewPoint->{1.3,-2.4,2}";}
289
290  str = str + closing + "];";
291
292  return(str);
293}
294example
295{ "EXAMPLE:"; echo =2;
296   // ---------  plane curves ------------
297   ring rr0 = 0,x1,dp; export rr0;
298
299   ideal I = x1^3 + x1, x1^2;
300   ideal J = x1^2,-x1+x1^3; "";
301   plot(J); "";
302   plot(I,"-2,2",J);
303
304   ring rr1 = 0,x,dp;  export rr1;
305   ideal I(1) = 2x2-1/2x3 +1, x-1;
306   ideal I(2) = x3-x ,x; "";
307   plot(I(1),I(2),"-2,4");
308
309 // Paste the output into a Mathematica notebook
310 // active evalutation of the cell with SHIFT RETURN
311
312 // Hit RETURN to continue
313 pause;
314   // --------- space curves --------------
315   I(1) = x3,-1/10x3+x2,x2;
316   I(2) = x2,-x2+2,-2x2-x+1; "";
317   plot(I(1),I(2));
318
319 // Paste the output into a Mathematica notebook
320 // active evalutation of the cell with SHIFT RETURN
321
322 // Hit RETURN to continue
323 pause;
324   // --------- implicit curves ------------
325   ring rr = 0,(x,y),ds;  export rr;
326   // implicit curves
327   ideal I(1) = y,-x2;
328   ideal I(2) = x2,-y2 +4; "";
329
330   //the lemniscate
331
332   ideal I(3) = x4+2x2y2 + y4, x2-y2;  "";
333
334   // a critical part
335   // adjust the plotregion properly to get a good picture
336
337   poly f = (x-y)*(x2+y);
338   plot(f); "";
339   ideal J = jacob(f);
340   J; "";
341   plot(J[1],J[2]);     // bad resolution
342   "";
343   plot(J[1],J[2],"-10,10","-10,10");  // get better picture in this way
344
345   plot(I(1),I(2),"-2.5,2.5"); "";
346
347   plot(I(3),"-2,2");
348
349 // Paste the output into a Mathematica notebook
350 // active evalutation of the cell with SHIFT RETURN
351
352 // Hit RETURN to continue
353 pause;
354   // ----------- surfaces -------------------
355   ideal J(1) = 3xy4 + 2xy2, x5y3 + x + y6,10x2;
356
357   // the Whitney-Umbrella
358
359   ideal J(2) = xy,y,x2; "";
360   plot(J(1),"-2,1","1,2");  "";
361   plot(J(2),"-1.5,1.5","-2,2");
362
363 // Paste the output into a Mathematica notebook
364 // active evalutation of the cell with SHIFT RETURN
365}
366///////////////////////////////////////////////////////////////////////////////
367
368
Note: See TracBrowser for help on using the repository browser.