1-- $Id: plot.cpkg5,v 1.2 2016/10/10 18:46:10 abbott Exp $
2Package $plot -- plotting functions
3
4Export ImplicitPlot;
5Export ImplicitPlotOn;
6Export PlotPoints;
7Export PlotPointsOn;
8
9Define About()
10  PrintLn  "    Author: J Abbott";
11  PrintLn  "    Date: 7 November 2008";
12EndDefine; -- About
13
14
15Define ImplicitPlot(F, Xrange, Yrange)
16  ImplicitPlotOn(F, Xrange, Yrange, "CoCoAPlot");
17EndDefine; -- ImplicitPlot
18
19Define PlotPoints(Points)
20  PlotPointsOn(Points, "CoCoAPlot");
21EndDefine; -- PlotPoints
22
23Define PlotPointsOn(Points, PlotFileName)
24  Foreach P In Points Do
25    If len(P)<>2 Then error("PointPlotOn: points must have 2 coordinates"); EndIf;
26  EndForeach;
27  PlotFile := OpenOFile(PlotFileName); //// ANNA2010-12-16, "w");
28  Print "Plotting points...";
29  Foreach P In Points Do
30    PrintLn DecimalStr(P[1]), "  ", DecimalStr(P[2]) On PlotFile;
31  EndForeach;
32  close(PlotFile);
33  PrintLn "100%";
34  PrintLn len(Points), " plotted points have been placed in the file "+PlotFileName;
35EndDefine; -- PlotPointsOn
36
37
38Define ImplicitPlotOn(F, Xrange, Yrange, PlotFileName)
39  If type(F) <> RINGELEM Then error("ImplicitPlot: first arg must be a polynomial"); EndIf;
40  P := owner(F);
41  VarInds := [I In 1..NumIndets(P) | deg(F,indet(P, I))> 0];
42  If len(VarInds) <> 2 Then
43    error("ImplicitPlot: polynomial must be genuinely bivariate");
44  EndIf;
45  X := indet(P,VarInds[1]);
46  Y := indet(P,VarInds[2]);
47  X0 := Xrange[1];
48  Xdelta := Xrange[2]-Xrange[1];
49  Y0 := Yrange[1];
50  Ydelta := Yrange[2]-Yrange[1];
51  If Xdelta <= 0 Or Ydelta <= 0 Then
52    error("ImplicitPlot: the X and Y ranges must both be of positive width");
53  EndIf;
54  GridSize := 256;
55  // Make the quantum values rather smaller than the deltas, otherwise
56  // intersections of graphs can appear to be in quite the wrong place!
57  // Perhaps a factor of 1/256 is overkill, but it probably does not cost Time.
58  Xquantum := Xdelta/(256*GridSize);
59  Yquantum := Ydelta/(256*GridSize);
60  PlotFile := OpenOFile(PlotFileName); //// ANNA2010-12-16, "w");
61  Print "Plotting points...";
62  NextPrint := 0.1;
63  Count := 0;
64  // Scan across columns...
65  For I := 0 To GridSize Do
66    If I > 2*NextPrint*GridSize Then
67      Print 100*NextPrint,"%...";
68      NextPrint := NextPrint+0.1;
69    EndIf;
70    Xval := X0+Xdelta*I/GridSize;
71    G := subst(F, X, Xval);
72    If G <> 0 And deg(G) > 0 Then
73      RRA := RealRootsApprox(G, Yquantum, Yrange);
74      Foreach Yval1 In RRA Do
75        PrintLn DecimalStr(Xval), "  ", DecimalStr(Yval1) On PlotFile;
76        Count := Count + 1;
77      EndForeach;
78    EndIf;
79  EndFor;
80  // Scan across rows...
81  For J := 0 To GridSize Do
82    If J+GridSize > 2*NextPrint*GridSize Then
83      Print 100*NextPrint,"%...";
84      NextPrint := NextPrint+0.1;
85    EndIf;
86    Yval := Y0+Ydelta*J/GridSize;
87    G := subst(F, Y, Yval);
88    If G <> 0 And deg(G) > 0 Then
89      RRA := RealRootsApprox(G, Xquantum, Xrange);
90      Foreach Xval1 In RRA Do
91        PrintLn DecimalStr(Xval1), "  ", DecimalStr(Yval) On PlotFile;
92        Count := Count + 1;
93      EndForeach;
94    EndIf;
95  EndFor;
96  close(PlotFile);
97  PrintLn "100%";
98  PrintLn Count, " plotted points have been placed in the file "+PlotFileName;
99EndDefine; -- ImplicitPlotOn
100
101
102EndPackage;
103