Curve that touches three points
Curve that touches three points is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Draw a curve that touches 3 points (1 starting point, 2 medium, 3 final point)
- Do not use functions of a library, implement the curve() function yourself
- coordinates:(x,y) starting point (10,10) medium point (100,200) final point (200,10)
Action!
INCLUDE "H6:REALMATH.ACT"
TYPE Point=[INT x,y]
PROC QuadraticCurve(Point POINTER p1,p2,p3 REAL POINTER a,b,c)
REAL x1,y1,x2,y2,x3,y3,x11,x22,x33,m,n,tmp1,tmp2,tmp3,tmp4,r1
IntToRealForNeg(-1,r1)
IntToRealForNeg(p1.x,x1) IntToRealForNeg(p1.y,y1)
IntToRealForNeg(p2.x,x2) IntToRealForNeg(p2.y,y2)
IntToRealForNeg(p3.x,x3) IntToRealForNeg(p3.y,y3)
RealMult(x1,x1,x11) ;x11=x1^2
RealMult(x2,x2,x22) ;x22=x2^2
RealMult(x3,x3,x33) ;x33=x3^2
RealSub(x1,x2,m) ;m=x1-x2
RealSub(x3,x2,n) ;n=x3-x2
RealMult(m,n,tmp1) ;tmp1=m*n
IF IsNegative(tmp1) THEN
RealMult(m,r1,tmp1)
RealAssign(tmp1,m) ;m=-m
FI
RealSub(y1,y2,tmp1) ;tmp1=y1-y2
RealMult(n,tmp1,tmp2) ;tmp2=n*(y1-y2)
RealSub(y3,y2,tmp1) ;tmp1=y3-y2
RealMult(m,tmp1,tmp3) ;tmp3=m*(y3-y2)
RealAdd(tmp2,tmp3,tmp1) ;tmp1=n*(y1-y2)+m*(y3-y2)
RealSub(x11,x22,tmp2) ;tmp2=x1^2-x2^2
RealMult(n,tmp2,tmp3) ;tmp3=n*(x1^2-x2^2)
RealSub(x33,x22,tmp2) ;tmp2=x3^2-x2^2
RealMult(m,tmp2,tmp4) ;tmp4=m*(x3^2-x2^2)
RealAdd(tmp3,tmp4,tmp2) ;tmp2=n*(x1^2-x2^2)+m*(x3^2-x2^2)
RealDiv(tmp1,tmp2,a) ;a=(n*(y1-y2)+m*(y3-y2)) / (n*(x1^2-x2^2)+m*(x3^2-x2^2))
RealSub(x33,x22,tmp1) ;tmp1=x3^2-x2^2
RealMult(tmp1,a,tmp2) ;tmp2=(x3^2-x2^2)*a
RealSub(y3,y2,tmp1) ;tmp1=y3-y2
RealSub(tmp1,tmp2,tmp3) ;tmp3=(y3-y2)-(x3^2-x2^2)*a
RealSub(x3,x2,tmp1) ;tmp1=x3-x2
RealDiv(tmp3,tmp1,b) ;b=((y3-y2)-(x3^2-x2^2)*a) / (x3-x2)
RealMult(a,x11,tmp1) ;tmp1=a*x1^2
RealMult(b,x1,tmp2) ;tmp2=b*x1
RealSub(y1,tmp1,tmp3) ;tmp3=y1-a*x1^2
RealSub(tmp3,tmp2,c) ;c=y1-a*x1^2-b*x1
RETURN
PROC DrawPoint(INT x,y)
Plot(x-2,y-2) DrawTo(x+2,y-2)
DrawTo(x+2,y+2) DrawTo(x-2,y+2)
DrawTo(x-2,y-2)
RETURN
INT FUNC Min(INT a,b)
IF a<b THEN RETURN (a) FI
RETURN (b)
INT FUNC Max(INT a,b)
IF a>b THEN RETURN (a) FI
RETURN (b)
INT FUNC CalcY(REAL POINTER a,b,c INT xi)
REAL xr,xr2,yr,tmp1,tmp2,tmp3
INT yi
IntToRealForNeg(xi,xr) ;xr=x
RealMult(xr,xr,xr2) ;xr2=x^2
RealMult(a,xr2,tmp1) ;tmp1=a*x^2
RealMult(b,xr,tmp2) ;tmp2=b*x
RealAdd(tmp1,tmp2,tmp3) ;tmp3=a*x^2+b*x
RealAdd(tmp3,c,yr) ;y3=a*x^2+b*x+c
yi=Round(yr)
RETURN (yi)
PROC DrawCurve(Point POINTER p1,p2,p3)
REAL a,b,c
INT xi,yi,minX,maxX
QuadraticCurve(p1,p2,p3,a,b,c)
DrawPoint(p1.x,p1.y)
DrawPoint(p2.x,p2.y)
DrawPoint(p3.x,p3.y)
minX=Min(p1.x,p2.x)
minX=Min(minX,p3.x)
maxX=Max(p1.x,p2.x)
maxX=Max(maxX,p3.x)
yi=CalcY(a,b,c,minX)
Plot(minX,yi)
FOR xi=minX TO maxX
DO
yi=CalcY(a,b,c,xi)
DrawTo(xi,yi)
OD
RETURN
PROC Main()
BYTE CH=$02FC,COLOR1=$02C5,COLOR2=$02C6
Point p1,p2,p3
Graphics(8+16)
Color=1
COLOR1=$0C
COLOR2=$02
p1.x=10 p1.y=10
p2.x=100 p2.y=180
p3.x=200 p3.y=10
DrawCurve(p1,p2,p3)
DO UNTIL CH#$FF OD
CH=$FF
RETURN- Output:
Screenshot from Atari 8-bit computer
Ada
Find center and radius of circle that touches the 3 points. Solve with simple linear algebra. In this case no division by zero.
with Ada.Text_Io;
with Ada.Numerics.Generic_Elementary_Functions;
procedure Three_Point_Circle is
type Real is new Float;
package Real_Math is
new Ada.Numerics.Generic_Elementary_Functions (Real);
package Real_Io is
new Ada.Text_Io.Float_Io (Real);
use Real_Io, Ada.Text_Io;
-- Point P1
X1 : constant Real := 10.0;
Y1 : constant Real := 10.0;
-- Point P2
X2 : constant Real := 100.0;
Y2 : constant Real := 200.0;
-- Point P3
X3 : constant Real := 200.0;
Y3 : constant Real := 10.0;
-- Point P4 - midpoint between P1 and P2
X4 : constant Real := (X1 + X2) / 2.0;
Y4 : constant Real := (Y1 + Y2) / 2.0;
S4 : constant Real := (Y2 - Y1) / (X2 - X1); -- Slope P1-P2
A4 : constant Real := -1.0 / S4; -- Slope P4-Center
-- Y4 = A4 * X4 + B4 <=> B4 = Y4 - A4 * X4
B4 : constant Real := Y4 - A4 * X4;
-- Point P5 - midpoint between P2 and P3
X5 : constant Real := (X2 + X3) / 2.0;
Y5 : constant Real := (Y2 + Y3) / 2.0;
S5 : constant Real := (Y3 - Y2) / (X3 - X2); -- Slope P2-P3
A5 : constant Real := -1.0 / S5; -- Slope P5-Center
-- Y5 = A5 * X5 + B5 <=> B5 = Y5 - A5 * X5
B5 : constant Real := Y5 - A5 * X5;
-- Find center
-- Y = A4 * X + B4 -- Line 1
-- Y = A5 * X + B5 -- Line 2
-- Solve for X:
-- A4 * X + B4 = A5 * X + B5
-- A4 * X - A5 * X = B5 - B4
-- X * (A4 - A5) = B5 - B4
-- X = (B5 - B4) / (A4 - A5)
Xc : constant Real := (B5 - B4) / (A4 - A5);
Yc : constant Real := A4 * Xc + B4;
-- Radius
R : constant Real := Real_Math.Sqrt ((X1 - Xc) ** 2 + (Y1 - Yc) ** 2);
begin
Real_Io.Default_Exp := 0;
Real_Io.Default_Aft := 1;
Put ("Center : "); Put ("("); Put (Xc); Put (", "); Put (Yc); Put (")"); New_Line;
Put ("Radius : "); Put (R); New_Line;
end Three_Point_Circle;
- Output:
Center : (105.0, 81.3) Radius : 118.8
ALGOL 68
As with the Ada and FreeBASIC translation of it, finds the centre and radius of a circle through the points.
BEGIN # draw a curve that passes through 3 points - translated from the Ada #
# sample - finds the centre and radius of a circle through the points #
REAL x1 = 10; # Point P1 #
REAL y1 = 10;
REAL x2 = 100 # Point P2 #;
REAL y2 = 200;
REAL x3 = 200; # Point P3 #
REAL y3 = 10;
REAL x4 = ( x1 + x2 ) / 2; # Point P4 - midpoint between P1 and P2 #
REAL y4 = ( y1 + y2 ) / 2;
REAL s4 = ( y2 - y1 ) / ( x2 - x1 ); # Slope P1-P2 #
REAL a4 = -1 / s4; # Slope P4-Centre #
# Y4 = A4 * X4 + B4 <=> B4 = Y4 - A4 * X4 #
REAL b4 = y4 - a4 * x4;
REAL x5 = ( x2 + x3 ) / 2; # Point P5 - midpoint between P2 and P3 #
REAL y5 = ( y2 + y3 ) / 2;
REAL s5 = ( y3 - y2 ) / ( x3 - x2 ); # Slope P2-P3 #
REAL a5 = -1.0 / s5; # Slope P5-Centre #
# Y5 = A5 * X5 + B5 <=> B5 = Y5 - A5 * X5 #
REAL b5 = y5 - a5 * x5;
# Find centre #
# Y = A4 * X + B4 -- Line 1 #
# Y = A5 * X + B5 -- Line 2 #
# Solve for X: #
# A4 * X + B4 = A5 * X + B5 #
# A4 * X - A5 * X = B5 - B4 #
# X * (A4 - A5) = B5 - B4 #
# X = (B5 - B4) / (A4 - A5) #
REAL xc = ( b5 - b4 ) / ( a4 - a5 );
REAL yc = a4 * xc + b4;
# Radius #
REAL r = sqrt( ( x1 - xc ) ^ 2 + ( y1 - yc ) ^ 2 );
print( ( "Centre : (", fixed( xc, -6, 1 ), ", ", fixed( yc, -6, 1 ), " )", newline ) );
print( ( "Radius : ", fixed( r, -6, 1 ), newline ) )
END- Output:
Centre : ( 105.0, 81.3 ) Radius : 118.8
AutoHotkey
QuadraticCurve(p1,p2,p3){ ; Y = aX^2 + bX + c
x1:=p1.1, y1:=p1.2, x2:=p2.1, y2:=p2.2, x3:=p3.1, y3:=p3.2
m:=x1-x2, n:=x3-x2, m:= ((m*n)<0?-1:1) * m
a:=(n*(y1-y2)+m*(y3-y2)) / (n*(x1**2 - x2**2) + m*(x3**2 - x2**2))
b:=((y3-y2) - (x3**2 - x2**2)*a) / (x3-x2)
c:=y1 - a*x1**2 - b*x1
return [a,b,c]
}
Examples:
P1 := [10,10], P2 := [100,200], P3 := [200,10]
v := QuadraticCurve(p1,p2,p3)
a := v.1, b:= v.2, c:= v.3
for i, X in [10,100,200]{
Y := a*X**2 + b*X + c ; Y = aX^2 + bX + c
res .= "[" x ", " y "]`n"
}
MsgBox % "Y = " a " X^2 " (b>0?"+":"") b " X " (c>0?"+":"") c " `n" res
- for plotting, use code from RosettaCode: Plot Coordinate Pairs
Outputs:
Y = -0.021111 X^2 +4.433333 X -32.222222 [10, 10.000000] [100, 200.000000] [200, 10.000000]
F#
This task uses Lagrange_Interpolation#F#
// Curve that touches three points. Nigel Galloway: September 13th., 2023
open Plotly.NET
let points=let a=LIF([10;100;200],[10;200;10]).Expression in [10.0..200.0]|>List.map(fun n->(n,(Evaluate.evaluate (Map.ofList ["x",n]) a).RealValue))
Chart.Point(points)|>Chart.show
- Output:
FreeBASIC
' Point P1
Dim As Double X1 = 10.0
Dim As Double Y1 = 10.0
' Point P2
Dim As Double X2 = 100.0
Dim As Double Y2 = 200.0
' Point P3
Dim As Double X3 = 200.0
Dim As Double Y3 = 10.0
' Point P4 - midpoint between P1 and P2
Dim As Double X4 = (X1 + X2) / 2.0
Dim As Double Y4 = (Y1 + Y2) / 2.0
Dim As Double S4 = (Y2 - Y1) / (X2 - X1) ' Slope P1-P2
Dim As Double A4 = -1.0 / S4 ' Slope P4-Center
' Y4 = A4 * X4 + B4 <=> B4 = Y4 - A4 * X4
Dim As Double B4 = Y4 - A4 * X4
' Point P5 - midpoint between P2 and P3
Dim As Double X5 = (X2 + X3) / 2.0
Dim As Double Y5 = (Y2 + Y3) / 2.0
Dim As Double S5 = (Y3 - Y2) / (X3 - X2) ' Slope P2-P3
Dim As Double A5 = -1.0 / S5 ' Slope P5-Center
' Y5 = A5 * X5 + B5 <=> B5 = Y5 - A5 * X5
Dim As Double B5 = Y5 - A5 * X5
' Find center
' Y = A4 * X + B4 ' Line 1
' Y = A5 * X + B5 ' Line 2
' Solve for X:
' A4 * X + B4 = A5 * X + B5
' A4 * X - A5 * X = B5 - B4
' X * (A4 - A5) = B5 - B4
' X = (B5 - B4) / (A4 - A5)
Dim As Double Xc = (B5 - B4) / (A4 - A5)
Dim As Double Yc = A4 * Xc + B4
' Radius
Dim As Double R = Sqr((X1 - Xc) ^ 2 + (Y1 - Yc) ^ 2)
Print Using "Center : (###.#, ###.#)"; Xc; Yc
Print Using "Radius : ###.#"; R
Sleep
Go
There are, of course, an infinity of curves which can be fitted to 3 points. The most obvious solution is to fit a quadratic curve (using Lagrange interpolation) and so that's what we do here.
As we're not allowed to use library functions to draw the curve, we instead divide the x-axis of the curve between successive points into equal segments and then join the resulting points with straight lines.
The resulting 'curve' is then saved to a .png file where it can be viewed with a utility such as EOG.
package main
import "github.com/fogleman/gg"
var p = [3]gg.Point{{10, 10}, {100, 200}, {200, 10}}
func lagrange(x float64) float64 {
return (x-p[1].X)*(x-p[2].X)/(p[0].X-p[1].X)/(p[0].X-p[2].X)*p[0].Y +
(x-p[0].X)*(x-p[2].X)/(p[1].X-p[0].X)/(p[1].X-p[2].X)*p[1].Y +
(x-p[0].X)*(x-p[1].X)/(p[2].X-p[0].X)/(p[2].X-p[1].X)*p[2].Y
}
func getPoints(n int) []gg.Point {
pts := make([]gg.Point, 2*n+1)
dx := (p[1].X - p[0].X) / float64(n)
for i := 0; i < n; i++ {
x := p[0].X + dx*float64(i)
pts[i] = gg.Point{x, lagrange(x)}
}
dx = (p[2].X - p[1].X) / float64(n)
for i := n; i < 2*n+1; i++ {
x := p[1].X + dx*float64(i-n)
pts[i] = gg.Point{x, lagrange(x)}
}
return pts
}
func main() {
const n = 50 // more than enough for this
dc := gg.NewContext(210, 210)
dc.SetRGB(1, 1, 1) // White background
dc.Clear()
for _, pt := range getPoints(n) {
dc.LineTo(pt.X, pt.Y)
}
dc.SetRGB(0, 0, 0) // Black curve
dc.SetLineWidth(1)
dc.Stroke()
dc.SavePNG("quadratic_curve.png")
}
J
NB. coordinates:(x,y) starting point (10,10) medium point (100,200) final point (200,10) X=: 10 100 200 Y=: 10 200 10 NB. matrix division computes polynomial coefficients NB. %. implements singular value decomposition NB. in other words, we can also get best fit polynomials of lower order. polynomial=: (Y %. (^/ ([: i. #)) X)&p. assert 10 200 10 -: polynomial X NB. test Filter=: (#~`)(`:6) Round=: adverb def '<.@:(1r2&+)&.:(%&m)' assert 100 120 -: 100 8 Round 123 NB. test, round 123 to nearest multiple of 100 and of 8 NB. libraries not permitted, character cell graphics are used. GRAPH=: 50 50 $ ' ' NB. is an array of spaces NB. place the axes GRAPH=: '-' [`(([:<0; i.@:#)@:])`]} GRAPH GRAPH=: '|' [`(([:<0;~i.@:#)@:])`]} GRAPH GRAPH=: '+' [`((<0;0)"_)`]} GRAPH NB. origin NB. clip the domain. EXES=: ((<:&(>./X) *. (<./X)&<:))Filter 5 * i. 200 WHYS=: polynomial EXES NB. draw the curve 1j1 #"1 |. 'X' [`((<"1 WHYS ;&>&:([: 1 Round %&5) EXES)"_)`]} GRAPH NB. were we to use a library: load'plot' 'title 3 point fit' plot (j. polynomial) i.201
Java
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Paths;
public final class CurveThatTouchesThreePoints {
public static void main(String[] args) throws IOException {
final double x1 = 10.0, y1 = 10.0; // point P1
final double x2 = 100.0, y2 = 200.0; // point P2
final double x3 = 200.0, y3 = 10.0; // point P3
final double x4 = ( x1 + x2 ) / 2.0; // x-coordinate of midpoint of line L1 between P1 and P2
final double y4 = ( y1 + y2 ) / 2.0; // y-coordinate of midpoint of line L1 between P1 and P2
final double slope4 = ( y2 - y1 ) / ( x2 - x1 ); // gradient of line LI
final double gradient4 = -1.0 / slope4; // gradient of line O1 orthogonal to L1
final double intercept4 = y4 - gradient4 * x4; // intercept of line O1 on y-axis
// Line O1 has equation: Y = gradient4 * X + intercept4
final double x5 = ( x2 + x3 ) / 2.0; // x-coordinate of midpoint of line L2 between P2 and P3
final double y5 = ( y2 + y3 ) / 2.0; // y-coordinate of midpoint of line L2 between P2 and P3
final double slope5 = ( y3 - y2 ) / ( x3 - x2 ); // gradient of line L2
final double gradient5 = -1.0 / slope5; // gradient of line O2 orthogonal to L2
final double intercept5 = y5 - gradient5 * x5; // intercept of line O2 on y-axis
// Line O2 has equation: Y = gradient5 * X + intercept5
// Solving the equations for lines O1 and O2
final double centreX = ( intercept5 - intercept4 ) / ( gradient4 - gradient5 );
final double centreY = gradient4 * centreX + intercept4;
final double radius = Math.sqrt(Math.pow(x1 - centreX, 2) + Math.pow(y1 - centreY, 2));
final int size = 300;
final int marginX = 50;
final int marginY = 70;
StringBuilder svgText = new StringBuilder();
svgText.append("<svg xmlns='http://www.w3.org/2000/svg'");
svgText.append(" width='" + size + "' height='" + size + "'>\n");
svgText.append("<rect width='100%' height='100%' fill='cyan'/>\n");
svgText.append("<circle r='" + radius + "' cx='" + ( centreX + marginX )
+ "' cy='" + ( size - centreY - marginY ) + "' stroke='red' stroke-width='3' fill='cyan'/>");
svgText.append("<line x1='0' y1='" + ( size - marginY ) + "' x2='"
+ size + "' y2='" + ( size - marginY ) + "' stroke='black' stroke-width='1'/>");
svgText.append("<line x1='" + marginX + "' y1='0' x2='" + marginX + "' y2='"
+ size + "' stroke='black' stroke-width='1'/>");
svgText.append("'\n</svg>\n");
String svg = svgText.toString();
Files.write(Paths.get("./CurveThreePointsJava.svg"), svg.getBytes());
}
}
- Output:
Media:CurveThreePointsJava.svg
Julia
To make things more specific, the example below finds the circle determined by the points. The curve is then the arc between the 3 points.
using Plots
struct Point; x::Float64; y::Float64; end
# Find a circle passing through the 3 points
const p1 = Point(10, 10)
const p2 = Point(100, 200)
const p3 = Point(200, 10)
const allp = [p1, p2, p3]
# set up problem matrix and solve.
# if (x - a)^2 + (y - b)^2 = r^2 then for some D, E, F, x^2 + y^2 + Dx + Ey + F = 0
# therefore Dx + Ey + F = -x^2 - y^2
v = zeros(Int, 3)
m = zeros(Int, 3, 3)
for row in 1:3
m[row, 1:3] .= [allp[row].x, allp[row].y, 1]
v[row] = -(allp[row].x)^2 - (allp[row].y)^2
end
q = (m \ v) # [-210.0, -162.632, 3526.32]
a, b, r = -q[1] / 2, -q[2] / 2, sqrt((q[1]^2/4) + q[2]^2/4 - q[3])
println("The circle with center at x = $a, y = $b and radius $r.")
x = a-r:0.25:a+r
y0 = sqrt.(r^2 .- (x .- a).^2)
plt = plot(x, y0 .+ b, color = :red)
plot!(x, b .- y0, color = :red)
scatter!([p.x for p in allp], [p.y for p in allp], markersize = r / 10)
- Output:
The circle with center at x = 105.0, y = 81.31578947368422 and radius 118.78948534384199.
Lambdatalk
We find a curve interpolating three points using a bezier algorithm. A bezier curve built on 3 points, p0, p1, p2 doesn't interpolate p1. We compute a new point q symetric of the middle of p0, p2 with respect to p1. The curve built on p0, q, p2 interpolates p0, p1, p2.
bezier interpolation of 3 points
p(t) = 1*p0(1-t)2 + 2*p1(1-t)t + 1*p2t2
{def interpol
{lambda {:p0 :p1 :p2 :t :u} // u =1-t
{+ {* 1 {A.get 0 :p0} :u :u}
{* 2 {A.get 0 :p1} :u :t}
{* 1 {A.get 0 :p2} :t :t}}
{+ {* 1 {A.get 1 :p0} :u :u}
{* 2 {A.get 1 :p1} :u :t}
{* 1 {A.get 1 :p2} :t :t}} }}
-> interpol
two useful functions
{def middle
{lambda {:p1 :p2} // compute the middle point of p1 and p2
{A.new
{/ {+ {A.get 0 :p1} {A.get 0 :p2}} 2}
{/ {+ {A.get 1 :p1} {A.get 1 :p2}} 2} }}}
-> middle
{def symetric // compute the symmetric point of p1 with respect to p2
{lambda {:p1 :p2}
{A.new
{- {* 2 {A.get 0 :p2}} {A.get 0 :p1} }
{- {* 2 {A.get 1 :p2}} {A.get 1 :p1} } }}}
-> symetric
computing the curve
{def curve
{lambda {:pol :n}
{S.map {{lambda {:p0 :p1 :p2 :n :i}
{interpol :p0 :p1 :p2 {/ :i :n} {- 1 {/ :i :n}}}
} {A.get 0 :pol} {A.get 1 :pol} {A.get 2 :pol} :n}
{S.serie -1 {+ :n 1}} }}}
-> curve
drawing a point
{def dot
{lambda {:pt}
{circle
{@ cx="{A.get 0 :pt}" cy="{A.get 1 :pt}" r="5"
stroke="#0ff" fill="transparent" stroke-width="2"}}}}
-> dot
defining points
{def P0 {A.new 150 180}} -> P0
{def P1 {A.new 300 250}} -> P1
{def P2 {A.new 150 330}} -> P2
{def P02 {middle {P0} {P2}}} -> P02
{def P20 {symetric {P02} {P1}}} -> P20
{def P10 {middle {P1} {P0}}} -> P10
{def P01 {symetric {P10} {P2}}} -> P01
{def P21 {middle {P2} {P1}}} -> P21
{def P12 {symetric {P21} {P0}}} -> P12
drawing points and curves
{svg {@ width="500" height="500" style="background:#444;"}
{polyline {@ points="{curve {A.new {P0} {P20} {P2}} 20}"
stroke="#f00" fill="transparent" stroke-width="4"}}
{polyline {@ points="{curve {A.new {P1} {P01} {P0}} 20}"
stroke="#0f0" fill="transparent" stroke-width="4"}}
{polyline {@ points="{curve {A.new {P2} {P12} {P1}} 20}"
stroke="#00f" fill="transparent" stroke-width="4"}}
{dot {P0}} {dot {P1}} {dot {P2}}
{dot {P02}} {dot {P20}}
{dot {P10}} {dot {P01}}
{dot {P21}} {dot {P12}}
}
See the result in http://lambdaway.free.fr/lambdawalks/?view=bezier_3
Mathematica /Wolfram Language
Built-in
pts = {{10, 10}, {100, 200}, {200, 10}};
cs = Circumsphere[pts]
Graphics[{PointSize[Large], Point[pts], cs}]
- Output:
Outputs the circle:
Sphere[{105, 1545/19}, (5 Sqrt[203762])/19]
and a graphical representation of the input points and the circle.
Alternate implementation
pts = {{10, 10}, {100, 200}, {200, 10}};
createCircle[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] :=
With[{a = Det[({{x1, y1, 1}, {x2, y2, 1}, {x3, y3, 1}})],
d = -Det[({{x1^2 + y1^2, y1, 1}, {x2^2 + y2^2, y2,
1}, {x3^2 + y3^2, y3, 1}})],
e = Det[({{x1^2 + y1^2, x1, 1}, {x2^2 + y2^2, x2, 1}, {x3^2 + y3^2,
x3, 1}})],
f = -Det[({{x1^2 + y1^2, x1, y1}, {x2^2 + y2^2, x2,
y2}, {x3^2 + y3^2, x3, y3}})]},
Circle[{-(d/(2 a)), -(e/(2 a))}, Sqrt[(d^2 + e^2)/(4 a^2) - f/a]]]
cs = createCircle[pts]
Graphics[{PointSize[Large], Point[pts], cs}]
- Output:
Outputs the circle:
Circle[{105, 1545/19}, (5 Sqrt[203762])/19]
and a graphical representation of the input points and the circle.
Nim
import imageman
type
FPoint = tuple[x, y: float]
FPoints3 = array[3, FPoint]
func lagrange(p: FPoints3; x: float): float =
(x-p[1].x) * (x-p[2].x) / (p[0].x-p[1].x) / (p[0].x-p[2].x) * p[0].y +
(x-p[0].x) * (x-p[2].x) / (p[1].x-p[0].x) / (p[1].x-p[2].x) * p[1].y +
(x-p[0].x) * (x-p[1].x) / (p[2].x-p[0].x) / (p[2].x-p[1].x) * p[2].y
func points(p: FPoints3; n: int): seq[Point] =
result.setLen(2 * n + 1)
var dx = (p[1].x - p[0].x) / float(n)
for i in 0..<n:
let x = p[0].x + dx * float(i)
result[i] = (x.toInt, p.lagrange(x).toInt)
dx = (p[2].x - p[1].x) / float(n)
for i in n..2*n:
let x = p[1].x + dx * float(i - n)
result[i] = (x.toInt, p.lagrange(x).toInt)
const N = 50
const P: FPoints3 =[(10.0, 10.0), (100.0, 200.0), (200.0, 10.0)]
var img = initImage[ColorRGBF](210, 210)
img.fill(ColorRGBF([float32 1, 1, 1])) # White background.
let color = ColorRGBF([float32 0, 0, 0]) # Black.
img.drawPolyline(closed = false, color, P.points(N))
img.savePNG("curve.png", compression = 9)
ooRexx
Version 1
/* REXX ***************************************************************
* Compute the polynome satisfying three given Points
**********************************************************************/
pl='(10,10) (100,200) (200,10)'
Do i=1 To 3
Parse Var pl '(' x.i ',' y.i ')' pl
s.i=x.i**2 x.i 1 y.i
End
Parse Value lingl() With a b c
If a<>0 Then
gl=a'*x**2'
Else
gl=''
If b>0 & gl<>'' Then b='+'||b
If b<>0 Then gl=gl||b'*x'
If c>0 & gl<>'' Then c='+'||c
If c<>0 Then gl=gl||c
Say 'y='gl
Say 'x / f(x) / y'
Do i=1 To 3
Say x.i '/' fun(x.i) '/' y.i
End
Exit
fun:
Parse Arg x
Return a*x**2+b*x+c
lingl: Procedure Expose s.
/************************************************* Version: 25.11.1996 *
* Lösung eines linearen Gleichungssystems
* 22.11.1996 PA neu
***********************************************************************/
Numeric Digits 12
Do i=1 to 3
l=s.i
Do j=1 By 1 While l<>''
Parse Var l a.1.i.j l
End
m=j-1
End
n=i-1
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.1.i.j,6,2)
End
Call dbg s
End
Do ie=1 To i-1
u=ie
v=ie+1
Do kk=ie To n
If a.u.kk.ie<>0 Then Leave
End
Select
When kk=ie Then Nop
When kk>n Then Call ex 'eine Katastrophe'
Otherwise Do
Do jj=1 To m
temp=a.u.ie.jj
a.u.ie.jj=a.u.kk.jj
a.u.kk.jj=temp
End
Do ip=1 To n
s=''
Do jp=1 To m
s=s format(a.u.ip.jp,12,2)
End
Call dbg s
End
End
End
Do i=1 To n
Do j=1 To m
If i<=ie Then
a.v.i.j=a.u.i.j
Else
a.v.i.j=a.u.i.j*a.u.ie.ie-a.u.i.ie*a.u.ie.j
End
End
Call dbg copies('-',70)
Do i=1 To n
Do j=1 To m
If a.v.i.j<>0 Then Leave
End
Select
When j=m Then Call ex 'Widersprü�chliches Gleichungssystem'
When j>m Then Call ex 'Gleichungen sind linear abhängig'
Otherwise Nop
End
End
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.v.i.j,12,2)
End
Call dbg s
End
End
n1=n+1
Do i=n To 1 By -1
i1=i+1
x.i=a.v.i.n1/a.v.i.i
sub=0
Do j=i+1 To n
sub=sub+a.v.i.j*x.j
End
x.i=x.i-sub/a.v.i.i
End
sol=''
Do i=1 To n
sol=sol x.i
End
Return sol
ex:
Say arg(1)
Exit
dbg: Return
- Output:
y=-0.021111111111*x**2+4.43333333333*x-32.2222222222 x / f(x) / y 10 / 10.0000000 / 10 100 / 200.000000 / 200 200 / 10.0000008 / 10
Version 2 using fraction arithmetic
/* REXX ***************************************************************
* Compute the polynome satisfying three given Points
**********************************************************************/
Numeric Digits 20
pl='(10,10) (100,200) (200,10)'
Do i=1 To 3
Parse Var pl '(' x.i ',' y.i ')' pl
s.i=x.i**2 x.i 1 y.i
End
abc=lingl()
a=abc[1]
b=abc[2]
c=abc[3]
If a~numerator<>0 Then
gl=a'*x**2'
Else
gl=''
If b~numerator<>0 Then gl=gl'+'||b'*x'
If c~numerator<>0 Then gl=gl'+'||c
o='y='gl
o=replr(o,'-(','+(-')
o=replr(o,'=-(','=(-')
o=replr(o,'=','=+')
Say o
Say 'x / f(x) / y'
Do i=1 To 3
Say x.i '/' fun(x.i) '/' y.i
End
Exit
fun:
Parse Arg x
Return a*x**2+b*x+c
lingl: Procedure Expose s.
/************************************************* Version: 25.11.1996 *
* Lösung eines linearen Gleichungssystems
* 22.11.1996 PA neu
***********************************************************************/
Numeric Digits 20
Do i=1 to 3
l=s.i
Do j=1 By 1 While l<>''
Parse Var l a.1.i.j l
fa.1.i.j=.fraction~new(a.1.i.j,1)
End
m=j-1
End
n=i-1
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.1.i.j,20)
End
Call dbg s
End
Do ie=1 To i-1
u=ie
v=ie+1
Do kk=ie To n
If a.u.kk.ie<>0 Then Leave
End
Select
When kk=ie Then Nop
When kk>n Then Call ex 'eine Katastrophe'
Otherwise Do
Do jj=1 To m
temp=a.u.ie.jj
a.u.ie.jj=a.u.kk.jj
a.u.kk.jj=temp
ftemp=fa.u.ie.jj
fa.u.ie.jj=fa.u.kk.jj
fa.u.kk.jj=ftemp
End
Do ip=1 To n
s=''
Do jp=1 To m
s=s format(a.u.ip.jp,20)
End
Call dbg s
End
End
End
Do i=1 To n
Do j=1 To m
If i<=ie Then Do
a.v.i.j=a.u.i.j
fa.v.i.j=fa.u.i.j
End
Else Do
a.v.i.j=a.u.i.j*a.u.ie.ie-a.u.i.ie*a.u.ie.j
fa.v.i.j=fa.u.i.j*fa.u.ie.ie-fa.u.i.ie*fa.u.ie.j
End
End
End
Call dbg copies('-',70)
Do i=1 To n
Do j=1 To m
If a.v.i.j<>0 Then Leave
End
Select
When j=m Then Call ex 'Widersprü�chliches Gleichungssystem'
When j>m Then Call ex 'Gleichungen sind linear abhängig'
Otherwise Nop
End
End
Do i=1 To n
s=''
Do j=1 To m
s=s format(a.v.i.j,20)
End
Call dbg s
End
End
n1=n+1
Do i=n To 1 By -1
x.i=a.v.i.n1/a.v.i.i
fx.i=fa.v.i.n1/fa.v.i.i
sub=0
fsub=.fraction~new(0,1)
Do j=i+1 To n
sub=sub+a.v.i.j*x.j
fsub=fsub+fa.v.i.j*fx.j
End
x.i=x.i-sub/a.v.i.i
fx.i=fx.i-fsub/fa.v.i.i
End
Return .array~of(fx.1,fx.2,fx.3)
ex:
Say arg(1)
Exit
dbg: Return
--REQUIRES fraction.cls
::class fraction public inherit stringlike orderable comparable
::method init /* initialize a fraction */
expose numerator denominator /* expose the state data */
Numeric Digits 20
Use Strict Arg numerator = 0, denominator = 1 /* access the two numbers */
numerator += 0 /* force rounding */
denominator += 0
anum=abs(numerator)
aden=abs(denominator)
x=gcd2(anum,aden)
anum=anum/x
aden=aden/x
If sign(denominator)<>sign(numerator) Then
numerator=-anum
Else
numerator=anum
denominator=aden
::method '[]' class /* create a new fraction */
forward message("NEW") /* just a synonym for NEW */
-- read-only attributes for numerator and denominator
::attribute numerator GET
::attribute denominator GET
::method '+' /* addition method */
expose numerator denominator /* access the state values */
Numeric Digits 20
Use Strict Arg adder = .nil /* get the operand */
if arg(1,'o') Then /* prefix plus operation? */
Return self /* don't do anything with this */
if adder~isa(.string) Then /* if just a simple number, */
adder = self~class~new(adder) /* convert to a fraction */
rnum=self~numerator*adder~denominator+,
self~denominator*adder~numerator
rdenom=self~denominator*adder~denominator
Return self~class~new(rnum,rdenom)
::method '-' /* subtraction method */
expose numerator denominator /* access the state values */
Numeric Digits 20
Use Strict Arg adder = .nil /* get the operand */
if arg(1,'o') Then do /* prefix minus operation? */
rdenom=self~denominator
rnum=-self~numerator
End
Else Do
if adder~isa(.string) Then /* if just a simple number, */
adder = self~class~new(adder) /* convert to a fraction */
rnum=self~numerator*adder~denominator-,
self~denominator*adder~numerator
rdenom=self~denominator*adder~denominator
End
Return self~class~new(rnum,rdenom)
::method '*' /* multiplication method */
expose numerator denominator /* access the state values */
Numeric Digits 20
Use Strict Arg adder = .nil /* get the operand */
if adder~isa(.string) Then /* if just a simple number, */
adder = self~class~new(adder) /* convert to a fraction */
rnum=self~numerator*adder~numerator
rdenom=self~denominator*adder~denominator
Return self~class~new(rnum,rdenom)
::method '/' /* division method */
expose numerator denominator /* access the state values */
Numeric Digits 20
Use Strict Arg adder = .nil /* get the operand */
if adder~isa(.string) Then /* if just a simple number, */
adder = self~class~new(adder) /* convert toa fraction */
rnum=self~numerator*adder~denominator
rdenom=self~denominator*adder~numerator
Return self~class~new(rnum,rdenom)
::method 'value' /* the fraction' numeric Value */
expose numerator denominator /* access the state values */
Return numerator/denominator
::method string /* format as a string value */
If self~denominator=1 Then
Return '('self~numerator')'
Else
Return '('self~numerator'/'self~denominator')' /* format as '(a,b)' */
::class "Stringlike" PUBLIC MIXINCLASS object
-- This unknown method forwards all method invocations to the object's string value,
-- effectively adding all of the string methods to the class
::method unknown UNGUARDED /* create an unknown method */
Use Arg msgname, args /* get the message and arguments */
/* just forward to the string val.*/
forward to(self~string) message(msgname) arguments(args)
::ROUTINE gcd2
/**********************************************************************
* Compute greatest common divider
**********************************************************************/
Numeric Digits 20
Parse Arg a,b
if b = 0 Then Return abs(a)
Return GCD2(b,a//b)
::ROUTINE replr
/* REXX ***************************************************************
* Replace,in s, occurrences of old by new and return the changed string
* ooRexx has the builtin function changestr
**********************************************************************/
Parse Arg s,new,old
Do i=1 To 2 Until p=0
p=pos(old,s)
If p>0 Then
s=left(s,p-1)||new||substr(s,p+length(old))
End
Return s
- Output:
y=-(19/900)*x**2+(133/30)*x-(290/9) x / f(x) / y 10 / (10) / 10 100 / (200) / 200
Version 3 computing the circumcircle (among many other things)
/* REXX ****************************************************************
* Triangle computes data about a given triangle
* The circumcircle is what we need here
***********************************************************************/
call triangle 10 10 200 10 100 200
Exit
triangle:
/***********************************************************************
* Triangle Computations
* 940810 PA new
* 220624 a mere 38 years later completed and anglisized
***********************************************************************/
Parse Arg ax ay bx by cx cy
If ax='?' Then Do
Say 'REXX Triangle ax ay bx by cx cy'
Say ' computes miscellaneous data about this triangle'
Exit
End
If ax='' Then Do
d='D 0 0 10 0 5 10'
Parse Var d . ax ay bx by cx cy .
End
Else
d='D' ax ay bx by cx cy .
Say ''
Say 'Triangle ABC:'
A='P' ax ay ; Say 'A' rep(A)
B='P' bx by ; Say 'B' rep(B)
C='P' cx cy ; Say 'C' rep(C)
areal=a(. ax ay bx by cx cy)
If areal<1e-3 Then
Call ex 'This isn''t a Triangle!! area='areal
Say ''
Say 'Triangle''s sides:'
al=dist(B,C) ; Say 'B-C a='round(al)
bl=dist(C,A) ; Say 'C-A b='round(bl)
cl=dist(A,B) ; Say 'A-B c='round(cl)
/* c**2=a**2+b**2-2*a*b*cos(gamma) */
cnvf=180/rxcalcpi() -- 57.2957796
alpha=rxCalcarccos((bl**2+cl**2-al**2)/(2*bl*cl),,'R')*cnvf
beta =rxCalcarccos((al**2+cl**2-bl**2)/(2*al*cl),,'R')*cnvf
gamma=rxCalcarccos((al**2+bl**2-cl**2)/(2*al*bl),,'R')*cnvf
Say ''
Say 'Triangle''s angles:'
Say 'alpha='round(alpha)
Say 'beta ='round(beta)
Say 'gamma='round(gamma)
Say 'sum ='round(alpha+beta+gamma)
Say ''
Say 'Angle-bisectors:'
wsa=ws(A,C,B); Say 'wsA' left(rep(wsA),20)
wsb=ws(B,A,C); Say 'wsB' left(rep(wsB),20)
wsc=ws(C,A,B); Say 'wsC' left(rep(wsC),20)
ha=normale(A,g(B,C))
Call dbg 'HA' rep(ha) ha
hb=normale(B,g(A,C))
Call dbg 'HB' rep(hb) hb
hc=normale(C,g(B,A))
Call dbg 'Hc' rep(hc) hc
HSP=sp(ha,hc)
If HSP='?' Then
HSP=sp(ha,hb)
Say ''
Say 'Orthocenter:' rep(HSP)
/***********************************************************************
* Perimeter and Area
***********************************************************************/
Say ''
Say 'Perimeter:' round(u(d))
Say 'Area: ' round(a(d))
/***********************************************************************
* Circumcircle
***********************************************************************/
U=sp(ss(A,B),ss(B,C))
Call dbg 'ss(A,B)='ss(A,B)
Call dbg 'ss(B,c)='ss(B,c)
Say ''
Say 'Center of circumcircle :' rep(U)
Say 'Radius :' round(dist(U,A))
/***********************************************************************
* Inscribed circle
***********************************************************************/
I=sp(wsa,wsb)
Say ''
Say 'Center of inscribed circle:' rep(I)
Say 'Radius :' round(rho(d))
/***********************************************************************
* Centroid
***********************************************************************/
Call dbg MP(B,C)
Call dbg MP(C,A)
sa=g(A,MP(B,C)); Call dbg 'sa='sa rep(sa)
sb=g(B,MP(C,A)); Call dbg 'sb='sb rep(sb)
S=sp(sa,sb)
Say ''
Say 'centroid:' rep(S)
/***********************************************************************
* Feuerbach Circle
***********************************************************************/
MAB='P' (ax+bx)/2 (ay+by)/2
MBC='P' (bx+cx)/2 (by+cy)/2
MCA='P' (cx+ax)/2 (cy+ay)/2
F=sp(ss(MAB,MBC),ss(MBC,MCA))
Say ''
Say 'Center of Feuerbach Circle:' rep(F)
Say 'Radius :' round(dist(F,MAB))
/***********************************************************************
* Euler's Line contains the following points:
* Centroid
* Center of circumcircle
* Orthocenter
* Center of Feuerbach Circle
***********************************************************************/
Call dbg 'Centroid..................' rep(S)
Call dbg 'Center of circumcircle....' rep(U)
Call dbg 'Orthocenter...............' rep(HSP)
Call dbg 'Center of Feuerbach Circle' rep(F)
Say ''
If abs(al-bl)<1.e-5 & abs(bl-cl)<1.e-5 Then
Say 'Equilateral Triangle - no Eulersche Gerade'
Else Do
Say 'Euler''s Line:'
su=rep(g(S,U)); Say 'S-U' su
sh=rep(g(S,HSP)); Say 'S-H' sh
sf=rep(g(S,F)); Say 'S-F' sf
uh=rep(g(U,HSP)); Say 'U-H' uh
End
Exit
round: Procedure
Numeric Digits 6
Parse Arg z
Return z+0
rep: Procedure Expose sigl
/***********************************************************************
* Show representation of a point or a line
***********************************************************************/
Parse Arg type a
Select
When type='P' Then Do
Parse Var a x y
res='('||round(x)||'/'||round(y)||')'
End
When type='g' Then Do
Parse Var a x1 y1 x2 y2
Select
When x1=x2 Then
res='x='||round(x1)
When y1=y2 Then
res='y='round(y1)
Otherwise Do
k=(y2-y1)/(x2-x1)
d=round(y1-k*x1)
Select
When d>0 Then d='+'d
When d=0 Then d=''
Otherwise Nop
End
If k=1 Then
res='y=x'd
Else
res='y='round(k)'*x'd
End
End
End
Otherwise Do
Say 'sigl='sigl
Say 'Unsupported type' type
res='???'
End
End
Return res
a: Procedure
/***********************************************************************
* Area (Heron's formula)
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
c=dist('P' ax ay,'P' bx by)
a=dist('P' bx by,'P' cx cy)
b=dist('P' cx cy,'P' ax ay)
s=(a+b+c)/2
res=rxCalcsqrt(s*(s-a)*(s-b)*(s-c))
Return res
rho: Procedure Expose ax ay bx by cx cy
/***********************************************************************
* Radius of inscribed circle
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
c=dist('P' ax ay,'P' bx by)
a=dist('P' bx by,'P' cx cy)
b=dist('P' cx cy,'P' ax ay)
s=(a+b+c)/2
res=rxCalcsqrt((s-a)*(s-b)*(s-c)/s)
Return res
u: Procedure
/***********************************************************************
* Perimeter
***********************************************************************/
Parse Arg . ax ay bx by cx cy .
Return dist('P' ax ay,'P' bx by)+,
dist('P' bx by,'P' cx cy)+,
dist('P' cx cy,'P' ax ay)
dist: Procedure
/***********************************************************************
* Distance between two points
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return rxCalcsqrt((x2-x1)**2+(y2-y1)**2)
g: Procedure
/***********************************************************************
* Intern representation of a line though two points
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return 'g' x1 y1 (x1+(x2-x1)) (y1+(y2-y1))
MP: Procedure
/***********************************************************************
* Center of a line segment
***********************************************************************/
Parse Arg . x1 y1 . , . x2 y2 .
Return 'P' ((x1+x2)/2) ((y2+y1)/2)
sp: Procedure
/***********************************************************************
* Intersection of two lines
***********************************************************************/
Parse Arg . xa ya xb yb . , . xc yc xd yd .
z=(xc-xa)*(yd-yc) - (yc-ya)*(xd-xc)
n=(xb-xa)*(yd-yc) - (yb-ya)*(xd-xc)
If n=0 Then Do
If z=0 Then
Call dbg 'lines are identical' z'/'n xa ya xb yb xc yc xd yd
Else
Call dbg 'lines are paralllel' z'/'n xa ya xb yb xc yc xd yd
Return '?'
End
Else Do
t=z/n
x=xa+(xb-xa)*t
y=ya+(yb-ya)*t
Call dbg x y
Return 'P' x y
End
euler: Procedure Expose S U HSP
/***********************************************************************
* Schwerpunkt, Umkreismittelpunkt, Höhenschnittpunkt
***********************************************************************/
Parse Arg . sx sy . ux uy . hx hy
Say 'Euler:' sx sy ux uy hx hy
eg=g(S,U); Say rep(eg)
eg2=g(S,HSP); Say rep(eg2)
eg3=g(U,HSP); Say rep(eg3)
Return
dist2:Procedure
/***********************************************************************
* Distance of a point C from a line AB
***********************************************************************/
Parse Arg ax ay bx by cx cy
Say 'A('ax'/'ay')' 'B('bx'/'by')' 'C('cx'/'cy')'
gx.1=ax
gx.2=bx-ax
gy.1=ay
gy.2=by-ay
Select
When gx.2=0 & gy.2=0 Then
Call ex 'g isn''t a line'
When gx.2=0 Then Do
xf=1
yf=0
c=-ax
End
When gy.2=0 Then Do
xf=0
yf=1
c=-ay
End
Otherwise Do
xf=1/gx.2
yf=-1/gy.2
c=-((ay/gy.2)+(ax/gx.2))
End
End
call dbg xf'*x+'yf'*y-'c'=0'
d=abs((xf*cx+yf*cy-c)/rxCalcsqrt(xf**2+yf**2))
Call dbg 'd='d
Return d
normale: Procedure
/***********************************************************************
* compute the line through point C that is normal to line A-B
***********************************************************************/
Parse Arg . ax ay . , . bx by cx cy .
vx=cx-bx
vy=cy-by
res='g' ax ay ax+vy ay-vx
Call dbg res
Return res
ss: Procedure
/***********************************************************************
* compute the perpendicular bisector of a line segment
***********************************************************************/
Parse Arg . ax ay . , . bx by .
Call dbg 'A('ax'/'ay')' 'B('bx'/'by')'
If ax=bx & ay=by Then
Call ex 'AB isn''t a line segment'
mx=(ax+bx)/2
my=(ay+by)/2
vx=bx-ax
vy=by-ay
Select
When vx=0 Then Parse Value 1 0 With sx sy
When vy=0 Then Parse Value 0 1 With sx sy
Otherwise Do
sx=vy
sy=-vx
End
End
Call dbg 'g' mx my (mx+sx) (my+sy)
Return 'g' mx my (mx+sx) (my+sy)
ws: Procedure
/***********************************************************************
* compute the angular symmetric of point A
***********************************************************************/
Parse Arg . ax ay . , . bx by . , . cx cy .
ebl=rxCalcsqrt((bx-ax)**2+(by-ay)**2)
ecl=rxCalcsqrt((cx-ax)**2+(cy-ay)**2)
--Say 'AB ' (bx-ax)/ebl (by-ay)/ebl
--Say 'AC ' (cx-ax)/ecl (cy-ay)/ecl
--Say 'AB+AC' ((bx-ax)/ebl+(cx-ax)/ecl) ((by-ay)/ebl+(cy-ay)/ecl)
res='g' ax ay ax+((bx-ax)/ebl+(cx-ax)/ecl)*10,
ay+((by-ay)/ebl+(cy-ay)/ecl)*10
Return res
dbg:
Return
Say arg(1)
ex:
Say arg(1)
Exit
::requires rxMath library
- Output:
Triangle ABC: A (10/10) B (200/10) C (100/200) Triangle's sides: B-C a=214.709 C-A b=210.238 A-B c=190 Triangle's angles: alpha=64.6538 beta =62.2415 gamma=53.1047 sum =180.000 Angle-bisectors: wsA y=0.632831*x+3.67169 wsB y=-0.603732*x+130.74 wsC y=-47.4947*x+4949.47 Orthocenter: (100.000/57.3684) Perimeter: 614.947 Area: 18050 Center of circumcircle : (105/81.3158) Radius : 118.789 Center of inscribed circle: (102.764/68.7042) Radius : 58.7042 centroid: (103.333/73.3333) Center of Feuerbach Circle: (102.500/69.3421) Radius : 59.3947 Euler's Line: S-U y=4.78947*x-421.579 S-H y=4.78947*x-421.579 S-F y=4.78948*x-421.579 U-H y=4.78947*x-421.579
Perl
Hilbert curve task code repeated here, with the addition that the 3 task-required points are marked. Satisfies the letter-of-the-law of task specification while (all in good fun) subverting the spirit of the thing.
use SVG;
use List::Util qw(max min);
use constant pi => 2 * atan2(1, 0);
# Compute the curve with a Lindemayer-system
%rules = (
A => '-BF+AFA+FB-',
B => '+AF-BFB-FA+'
);
$hilbert = 'A';
$hilbert =~ s/([AB])/$rules{$1}/eg for 1..6;
# Draw the curve in SVG
($x, $y) = (0, 0);
$theta = pi/2;
$r = 5;
for (split //, $hilbert) {
if (/F/) {
push @X, sprintf "%.0f", $x;
push @Y, sprintf "%.0f", $y;
$x += $r * cos($theta);
$y += $r * sin($theta);
}
elsif (/\+/) { $theta += pi/2; }
elsif (/\-/) { $theta -= pi/2; }
}
$max = max(@X,@Y);
$xt = -min(@X)+10;
$yt = -min(@Y)+10;
$svg = SVG->new(width=>$max+20, height=>$max+20);
$points = $svg->get_path(x=>\@X, y=>\@Y, -type=>'polyline');
$svg->rect(width=>"100%", height=>"100%", style=>{'fill'=>'black'});
$svg->polyline(%$points, style=>{'stroke'=>'orange', 'stroke-width'=>1}, transform=>"translate($xt,$yt)");
my $task = $svg->group( id => 'task-points', style => { stroke => 'red', fill => 'red' },);
$task->circle( cx => 10, cy => 10, r => 1, id => 'point1' );
$task->circle( cx => 100, cy => 200, r => 1, id => 'point2' );
$task->circle( cx => 200, cy => 10, r => 1, id => 'point3' );
open $fh, '>', 'curve-3-points.svg';
print $fh $svg->xmlify(-namespace=>'svg');
close $fh;
Phix
You can run this online here.
with javascript_semantics include pGUI.e Ihandle dlg, canvas cdCanvas cddbuffer, cdcanvas enum X, Y constant p = {{10,10},{100,200},{200,10}} function lagrange(atom x) return (x - p[2][X])*(x - p[3][X])/(p[1][X] - p[2][X])/(p[1][X] - p[3][X])*p[1][Y] + (x - p[1][X])*(x - p[3][X])/(p[2][X] - p[1][X])/(p[2][X] - p[3][X])*p[2][Y] + (x - p[1][X])*(x - p[2][X])/(p[3][X] - p[1][X])/(p[3][X] - p[2][X])*p[3][Y] end function function getPoints(integer n) sequence pts = {} atom {dx,pt,cnt} := {(p[2][X] - p[1][X])/n, p[1][X], n} for j=1 to 2 do for i=0 to cnt do atom x := pt + dx*i; pts = append(pts,{x,lagrange(x)}); end for {dx,pt,cnt} = {(p[3][X] - p[2][X])/n, p[2][X], n+1}; end for return pts end function procedure draw_cross(sequence xy) integer {x,y} = xy cdCanvasLine(cddbuffer, x-3, y, x+3, y) cdCanvasLine(cddbuffer, x, y-3, x, y+3) end procedure function redraw_cb(Ihandle /*ih*/) cdCanvasActivate(cddbuffer) cdCanvasSetForeground(cddbuffer, CD_BLUE) cdCanvasBegin(cddbuffer,CD_OPEN_LINES) atom {x,y} = {p[1][X], p[1][Y]}; -- curve starting point cdCanvasVertex(cddbuffer, x, y) sequence pts = getPoints(50) for i=1 to length(pts) do {x,y} = pts[i] cdCanvasVertex(cddbuffer, x, y) end for cdCanvasEnd(cddbuffer) cdCanvasSetForeground(cddbuffer, CD_RED) for i=1 to length(p) do draw_cross(p[i]) end for cdCanvasFlush(cddbuffer) return IUP_DEFAULT end function function map_cb(Ihandle ih) cdcanvas = cdCreateCanvas(CD_IUP, ih) cddbuffer = cdCreateCanvas(CD_DBUFFER, cdcanvas) cdCanvasSetBackground(cddbuffer, CD_WHITE) return IUP_DEFAULT end function procedure main() IupOpen() canvas = IupCanvas(NULL) IupSetAttribute(canvas, "RASTERSIZE", "220x220") IupSetCallback(canvas, "MAP_CB", Icallback("map_cb")) dlg = IupDialog(canvas,"DIALOGFRAME=YES") IupSetAttribute(dlg, "TITLE", "Quadratic curve") IupSetCallback(canvas, "ACTION", Icallback("redraw_cb")) IupMap(dlg) IupSetAttribute(canvas, "RASTERSIZE", NULL) IupShowXY(dlg,IUP_CENTER,IUP_CENTER) if platform()!=JS then IupMainLoop() IupClose() end if end procedure main()
PHP
Finds a circle that passes through those three points. Adapted from AutoHotkey.
function circle_3_points($x1, $y1, $x2, $y2, $x3, $y3) {
$x4 = ($x1 + $x2) / 2.0;
$y4 = ($y1 + $y2) / 2.0;
$s4 = ($y2 - $y1) / ($x2 - $x1);
$a4 = -1.0 / $s4;
$b4 = $y4 - $a4 * $x4;
$x5 = ($x2 + $x3) / 2.0;
$y5 = ($y2 + $y3) / 2.0;
$s5 = ($y3 - $y2) / ($x3 - $x2);
$a5 = -1.0 / $s5;
$b5 = $y5 - $a5 * $x5;
$xc = ($b5 - $b4) / ($a4 - $a5);
$yc = $a4 * $xc + $b4;
$r = √(²($x1 - $xc) + ²($y1 - $yc));
return [$xc, $yc, $r];
}
Raku
(formerly Perl 6)
Kind of bogus. There are an infinite number of curves that pass through those three points. I'll assume a quadratic curve. Lots of bits and pieces borrowed from other tasks to avoid relying on library functions.
Saved as a png for wide viewing support. Note that png coordinate systems have 0,0 in the upper left corner.
use Image::PNG::Portable;
# the points of interest
my @points = (10,10), (100,200), (200,10);
# Solve for a quadratic line that passes through those points
my (\a, \b, \c) = (rref ([.[0]², .[0], 1, .[1]] for @points) )[*;*-1];
# Evaluate quadratic equation
sub f (\x) { a×x² + b×x + c }
my ($w, $h) = 500, 500; # image size
my $scale = 2; # scaling factor
my $png = Image::PNG::Portable.new: :width($w), :height($h);
my ($lastx, $lasty) = 8, f(8).round;
(9 .. 202).map: -> $x {
my $f = f($x).round;
line($lastx, $lasty, $x, $f, $png, [0,255,127]);
($lastx, $lasty) = $x, $f;
}
# Highlight the defining points
dot( | $_, $(255,0,0), $png, 2) for @points;
$png.write: 'Curve-3-points-perl6.png';
# Assorted helper routines
sub rref (@m) {
return unless @m;
my ($lead, $rows, $cols) = 0, @m, @m[0];
for ^$rows -> $r {
$lead < $cols or return @m;
my $i = $r;
until @m[$i;$lead] {
++$i == $rows or next;
$i = $r;
++$lead == $cols and return @m;
}
@m[$i, $r] = @m[$r, $i] if $r != $i;
@m[$r] »/=» $ = @m[$r;$lead];
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
}
++$lead;
}
@m
}
sub line($x0 is copy, $y0 is copy, $x1 is copy, $y1 is copy, $png, @rgb) {
my $steep = abs($y1 - $y0) > abs($x1 - $x0);
($x0,$y0,$x1,$y1) »×=» $scale;
if $steep {
($x0, $y0) = ($y0, $x0);
($x1, $y1) = ($y1, $x1);
}
if $x0 > $x1 {
($x0, $x1) = ($x1, $x0);
($y0, $y1) = ($y1, $y0);
}
my $Δx = $x1 - $x0;
my $Δy = abs($y1 - $y0);
my $error = 0;
my $Δerror = $Δy / $Δx;
my $y-step = $y0 < $y1 ?? 1 !! -1;
my $y = $y0;
next if $y < 0;
for $x0 .. $x1 -> $x {
next if $x < 0;
if $steep {
$png.set($y, $x, |@rgb);
} else {
$png.set($x, $y, |@rgb);
}
$error += $Δerror;
if $error ≥ 0.5 {
$y += $y-step;
$error -= 1.0;
}
}
}
sub dot ($X is copy, $Y is copy, @rgb, $png, $radius = 3) {
($X, $Y) »×=» $scale;
for ($X X+ -$radius .. $radius) X ($Y X+ -$radius .. $radius) -> ($x, $y) {
$png.set($x, $y, |@rgb) if ( $X - $x + ($Y - $y) × i ).abs <= $radius;
}
}
See Curve-3-points-perl6.png (offsite .png image)
RPL
HP-49+ RPL has a built-in function named LAGRANGE that returns the curve equation from the 3 points, but let's consider it as belonging to a library.
« → a b c
« { 0 0 0 }
c a - C→R SWAP / c b - RE /
b a - C→R SWAP / c b - RE / -
1 SWAP PUT
b a - C→R SWAP / OVER 1 GET b a + RE * -
2 SWAP PUT
a IM OVER 1 GET a RE SQ * - OVER 2 GET a RE * -
3 SWAP PUT
{ 'X^2' 'X' 1 } * ∑LIST
» » 'PAREQ' STO
« # 131d # 64d PDIM 0 210 DUP2 XRNG YRNG
FUNCTION
(10,10) (100,200) (200,10)
3 DUPN PAREQ STEQ
ERASE DRAW
1 3 START PIXOFF NEXT @ switch off the 3 points on the curve
{ } PVIEW
{ PPAR EQ } PURGE
» 'TASK' STO
Wren
import "graphics" for Canvas, Color, Point
import "dome" for Window
class Game {
static init() {
Window.title = "Quadratic curve"
var width = 210
var height = 210
Window.resize(width, height)
Canvas.resize(width, height)
Canvas.cls(Color.white)
var n = 50
var p = [Point.new(10, 10), Point.new(100, 200), Point.new(200, 10)]
var col = Color.black // black curve
quadratic(n, p, col)
}
static update() {}
static draw(alpha) {}
static lagrange(p, x) {
return (x-p[1].x)*(x-p[2].x)/(p[0].x-p[1].x)/(p[0].x-p[2].x)*p[0].y +
(x-p[0].x)*(x-p[2].x)/(p[1].x-p[0].x)/(p[1].x-p[2].x)*p[1].y +
(x-p[0].x)*(x-p[1].x)/(p[2].x-p[0].x)/(p[2].x-p[1].x)*p[2].y
}
static quadratic(n, p, col) {
var pts = List.filled(2*n+1, null)
var dx = (p[1].x - p[0].x) / n
for (i in 0...n) {
var x = p[0].x + dx*i
pts[i] = Point.new(x, lagrange(p, x))
}
dx = (p[2].x - p[1].x) / n
for (i in n...2*n+1) {
var x = p[1].x + dx*(i-n)
pts[i] = Point.new(x, lagrange(p, x))
}
var prev = pts[0]
for (pt in pts.skip(1)) {
Canvas.line(prev.x, prev.y, pt.x, pt.y, col)
prev = pt
}
}
}
XPL0
def X0=10., Y0=10., X1=100., Y1=200., X2=200., Y2=10.;
func real Lagrange(X);
real X;
return (X-X1) * (X-X2) / (X0-X1) / (X0-X2) * Y0 +
(X-X0) * (X-X2) / (X1-X0) / (X1-X2) * Y1 +
(X-X0) * (X-X1) / (X2-X0) / (X2-X1) * Y2;
def Offset = 205;
real X;
[SetVid($13); \VGA 320x200x8 graphics
X:= X0;
Move(fix(X), Offset-fix(Y0));
repeat X:= X + 1.;
Line(fix(X), Offset-fix(Lagrange(X)), 3\cyan\);
until X >= X2;
Point(fix(X0), Offset-fix(Y0), 14\yellow\);
Point(fix(X1), Offset-fix(Y1), 14);
Point(fix(X2), Offset-fix(Y2), 14);
]zkl
Uses Image Magick and the PPM class from http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#zkl
const X=0, Y=1; // p.X == p[X]
var p=L(L(10.0, 10.0), L(100.0, 200.0), L(200.0, 10.0)); // (x,y)
fcn lagrange(x){ // float-->float
(x - p[1][X])*(x - p[2][X])/(p[0][X] - p[1][X])/(p[0][X] - p[2][X])*p[0][Y] +
(x - p[0][X])*(x - p[2][X])/(p[1][X] - p[0][X])/(p[1][X] - p[2][X])*p[1][Y] +
(x - p[0][X])*(x - p[1][X])/(p[2][X] - p[0][X])/(p[2][X] - p[1][X])*p[2][Y]
}
fcn getPoints(n){ // int-->( (x,y) ..)
pts:=List.createLong(2*n+1);
dx,pt,cnt := (p[1][X] - p[0][X])/n, p[0][X], n;
do(2){
foreach i in (cnt){
x:=pt + dx*i;
pts.append(L(x,lagrange(x)));
}
dx,pt,cnt = (p[2][X] - p[1][X])/n, p[1][X], n+1;
}
pts
}
fcn main{
var [const] n=50; // more than enough for this
img,color := PPM(210,210,0xffffff), 0; // white background, black curve
foreach x,y in (p){ img.cross(x.toInt(),y.toInt(), 0xff0000) } // mark 3 pts
a,b := p[0][X].toInt(), p[0][Y].toInt(); // curve starting point
foreach x,y in (getPoints(n)){
x,y = x.toInt(),y.toInt();
img.line(a,b, x,y, color); // can only deal with ints
a,b = x,y;
}
img.writeJPGFile("quadraticCurve.zkl.jpg");
}();- Output:
Image at quadratic curve
