using Plots
using Printf
using SymPy
##################
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, lw=0.5)
plot!([x1, x2], [y1, y2], color=color; linestyle=linestyle, lw=lw)
end;
##################
function abline(x0, y0, slope, minx, maxx, color=:black; lw=0.5)
f(x) = slope * (x - x0) + y0
segment(minx, f(minx), maxx, f(maxx), color; lw=lw)
end;
##################
function intersectionXY(x1, y1, x2, y2, x3, y3, x4, y4)
denominator = (x1*y3 - x1*y4 - x2*y3 + x2*y4 - x3*y1 + x3*y2 + x4*y1 - x4*y2)
X = (x1*x3*y2 - x1*x3*y4 - x1*x4*y2 + x1*x4*y3 - x2*x3*y1 + x2*x3*y4 + x2*x4*y1 - x2*x4*y3) / denominator
Y = (x1*y2*y3 - x1*y2*y4 - x2*y1*y3 + x2*y1*y4 - x3*y1*y4 + x3*y2*y4 + x4*y1*y3 - x4*y2*y3) / denominator
(X, Y)
end
##################
function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2 # 二乗距離を返す！！
end;
function dist(x1, y1, x2, y2, x0, y0; point=false)
# @syms dx, dy, line_length, vx, vy, projection_length, nearest_point_x, nearest_point_y
# 直線の方向ベクトル
dx = x2 - x1
dy = y2 - y1
# 直線の長さ
line_length = sqrt(dx^2 + dy^2)
# 直線上の点までのベクトル
vx = x0 - x1
vy = y0 - y1
# 直線上の点までの射影ベクトルの長さ
projection_length = (vx * dx + vy * dy) / line_length
# 直線上の最近点の座標
nearest_point_x = x1 + (projection_length * dx) / line_length
nearest_point_y = y1 + (projection_length * dy) / line_length
if point
# 直線上の最近点の座標を返す
return (nearest_point_x, nearest_point_y)
else
# 点 (x0, y0) から直線までの二乗距離を返す
return (x0 - nearest_point_x)^2 + (y0 - nearest_point_y)^2
end
end;
function dist2(x1, y1, x2, y2, x0, y0, r)
@syms dummy_temp
return numerator(apart(dist(x1, y1, x2, y2, x0, y0) - r^2, dummy_temp)) |> simplify# |> simplify
end;
#### 垂線の脚の座標
foot(x1, y1, x2, y2, x0, y0) = dist(x1, y1, x2, y2, x0, y0; point=true)
##################
#### 射影（回転）
function transform(x, y; deg=0, dx=0, dy=0)
(x, y) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [x, y]
return (x .+ dx, y .+ dy)
end;
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, by=0.5, n=0, lw=0.5)
n != 0 && (by = (endangle - beginangle) / n)
θ = beginangle:by:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=lw)
end;
#### 回転角度を指定して複数個描く
function rotate(ox, oy, r, color=:red; angle=120, beginangle=0, endangle=360, by=0.5, n=0)
for deg in 0:angle:360-1
(ox2, oy2) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [ox; oy]
circle(ox2, oy2, r, color; beginangle, endangle, by, n)
beginangle += angle
endangle += angle
end
end;
function rotatef(ox, oy, r, color=:red; angle=120, beginangle=0, endangle=360, by=0.5, n=0)
for deg in 0:angle:360-1
(ox2, oy2) = [cosd(deg) sind(deg); -sind(deg) cosd(deg)] * [ox; oy]
circlef(ox2, oy2, r, color; beginangle, endangle, by, n)
end
end;
##### 2 個描く
function circle2(x, y, r, color=:red)
circle(x, y, r, color)
circle(-x, y, r, color)
end;
##### 2 個描く
function circle2f(x, y, r, color=:red)
circlef(x, y, r, color)
circlef(-x, y, r, color)
end;
##### 2 個描く --2
function circle22(x, y, r, color=:red)
circle(x, y, r, color)
circle(x, -y, r, color)
end;
##### 2 個描く --2
function circle22f(x, y, r, color=:red)
circlef(x, y, r, color)
circlef(x, -y, r, color)
end;
##### 4 個描く
function circle4(x, y, r, color=:red)
circle(x, y, r, color)
circle(x, -y, r, color)
circle(-x, y, r, color)
circle(-x, -y, r, color)
end;
##### 4 個描く
function circle4f(x, y, r, color=:red)
circlef(x, y, r, color)
circlef(x, -y, r, color)
circlef(-x, y, r, color)
circlef(-x, -y, r, color)
end;
##### 4 個描く --2
function circle42(x, y, r, color=:red)
circle(x, y, r, color)
circle(x, -y, r, color)
circle(y, -x, r, color)
circle(-y, -x, r, color)
end;
##### 4 個描く --2
function circle42f(x, y, r, color=:red)
circlef(x, y, r, color)
circlef(x, -y, r, color)
circlef(y, -x, r, color)
circlef(-y, -x, r, color)
end;
# 塗りつぶしバージョン
function circlef(ox, oy, r, color=:red; beginangle=0, endangle=360, by=0.5, n=0)
n != 0 && (by = (endangle - beginangle) / n)
θ = beginangle:by:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end;
# アノテーション
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true, delta=0, deltax=0)
mark && scatter!([x], [y], color=color, markerstrokewidth=0, markersize=3)
annotate!(x + deltax, y + delta, text(string, 10, position, color, vertical))
end;
##################
function dimension_line(x1, y1, x2, y2, str="", color=:black, position=:center, vertical=:vcenter; linestyle=:solid, lw=0.5, deltax=0, delta=0)
plot!([x1, x2], [y1, y2], color=color; arrow=:arrow, linestyle=linestyle, lw=lw)
plot!([x2, x1], [y2, y1], color=color; arrow=:arrow, linestyle=linestyle, lw=lw)
annotate!((x1 + x2)/2 + deltax, (y1 + y2)/2 + delta, text(str, position, color, vertical, 10))
end;
# 矩形
function rect(x1, y1, x2, y2, color=:pink; fill=false)
if fill == false
plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], color=color, linewidth=0.5)
else
plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end
end;
# 楕円
function ellipse(ox, oy, ra, rb; φ=0, beginangle=0, endangle=360,
color=:black, lty=:solid, lwd=0.5, fcolor="")
"""
(ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。
fcolor を指定すると塗りつぶし。
"""
θ = beginangle:0.1:endangle
if φ == 0
if fcolor == ""
plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
linecolor=color, linestyle=lty, linewidth=lwd)
else
plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
linecolor=color, linestyle=lty, linewidth=lwd,
seriestype=:shape, fillcolor=fcolor)
end
else
x = ra .* cosd.(θ)
y = rb .* sind.(θ)
cosine = cosd(φ)
sine = sind(φ)
if fcolor == ""
plot!(cosine .* x .- sine .* y .+ ox,
sine .* x .+ cosine .* y .+ oy,
linecolor=color, linestyle=lty, linewidth=lwd)
else
plot!(cosine .* x .- sine .* y .+ ox,
sine .* x .+ cosine .* y .+ oy,
linecolor=color, linestyle=lty, linewidth=lwd,
seriestype=:shape, fillcolor=fcolor)
end
end
end;
#=
引数
楕円の長径 a，短径 b，接線の傾き（符号に注意）
戻り値
接点の座標 (x, y)，楕円に外接する円の半径 r
func(a, b, tanθ) = (-a^2*tanθ/sqrt(a^2*tanθ^2 + b^2), b^2/sqrt(a^2*tanθ^2 + b^2), -b + sqrt(a^2*tanθ^2 + b^2));
=#;
##### 多角形の面積
function area(xy)
x = xy[:, 1]
sum((vcat(x[2:end], x[1]) - vcat(x[end], x[1:end-1])) .* xy[:, 2]) / 2
end;
##### 二直線の交点座標
function intersection(x1, y1, x2, y2, x3, y3, x4, y4)
p1, p2, p3, p4 = sympy.Point(x1, y1), sympy.Point(x2, y2), sympy.Point(x3, y3), sympy.Point(x4, y4)
l1, l2 = sympy.Line(p1, p2), sympy.Line(p3, p4)
a = l1.intersection(l2)[1]
(a.x, a.y)
end;
##### 楕円に内接する円の中心座標
function getd(a, b, r; x0 = 0, y0 = 0)
d = sqrt((a^2 - b^2)*(b^2 - r^2))/b
x = a^2*d/(a^2 - b^2)
y = sqrt(r^2 - (x - d)^2)
return (d, x, y)
end;
##### 楕円に内接する円の半径
function getr(a, b, d; x0 = 0, y0 = 0)
r = b*sqrt((a^2 - b^2 - d^2)/((a^2 - b^2)))
x = a^2*d/(a^2 - b^2)
y = sqrt(r^2 - (x - d)^2)
return (r, x, y)
end;
##### 曲率円の半径
geth(a, b) = b^2/a;
##### 外心を計算する関数
function circumcenter(x1, y1, x2, y2, x3, y3)
D = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
x_o = ((x1^2 + y1^2) * (y2 - y3) + (x2^2 + y2^2) * (y3 - y1) + (x3^2 + y3^2) * (y1 - y2)) / D
y_o = ((x1^2 + y1^2) * (x3 - x2) + (x2^2 + y2^2) * (x1 - x3) + (x3^2 + y3^2) * (x2 - x1)) / D
r = sqrt((x1 - x_o)^2 + (y1 - y_o)^2)
return r, x_o, y_o
end;
##### 内心を計算する関数
function incenter(x1, y1, x2, y2, x3, y3)
a = sqrt((x2 - x3)^2 + (y2 - y3)^2)
b = sqrt((x1 - x3)^2 + (y1 - y3)^2)
c = sqrt((x1 - x2)^2 + (y1 - y2)^2)
x_i = (a * x1 + b * x2 + c * x3) / (a + b + c)
y_i = (a * y1 + b * y2 + c * y3) / (a + b + c)
r = sqrt(distance(x1, y1, x2, y2, x_i, y_i))
return r, x_i, y_i
end;
