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