using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve( (vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve( (vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
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;
#### 回転角度を指定して複数個描く
#### 回転の中心座標 (center_x, center_y) を指定できるようにした 25/07/12
function rotate(ox, oy, r, color=:red; center_x=0, center_y=0, angle=120, beginangle=0, endangle=360, lw=0.5, by=0.5, n=0)
ox -= center_x
oy -= center_y
for deg in 0:angle:360-1
(ox2, oy2) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [ox; oy]
circle(ox2 + center_x, oy2 + center_y, r, color; beginangle, endangle, by, n, lw=0.5)
beginangle += angle
endangle += angle
end
end;
function rotatef(ox, oy, r, color=:red; center_x=0, center_y=0, angle=120, beginangle=0, endangle=360, lw=0.5, by=0.5, n=0)
ox -= center_x
oy -= center_y
for deg in 0:angle:360-1
(ox2, oy2) = [cosd(deg) sind(deg); -sind(deg) cosd(deg)] * [ox; oy]
circlef(ox2 + center_x, oy2 + center_y, 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; fontsize=10, mark=true, delta=0, deltax=0)
mark && scatter!([x], [y], color=color, markerstrokewidth=0, markersize=3)
annotate!(x + deltax, y + delta, text(string, fontsize, 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 dimension_line(x1, y1, x2, y2, str="", color=:black, position=:center, vertical=:vcenter; dx=0, dy=0, length=0.2, α=π/12, lw=0.5, deltax=0, delta=0)
# 寸法線（中心の直線部分）
plot!(dx .+ [x1, x2], dy .+ [y1, y2], color=color, lw=lw, label="")
# 矢印の向きを求める
θ = atan(y2 - y1, x2 - x1) # 線分の傾き
# 矢印の端点を計算
arrow1_x = x1 + length * cos(θ - α)
arrow1_y = y1 + length * sin(θ - α)
arrow2_x = x1 + length * cos(θ + α)
arrow2_y = y1 + length * sin(θ + α)
arrow3_x = x2 + length * cos(θ + π - α)
arrow3_y = y2 + length * sin(θ + π - α)
arrow4_x = x2 + length * cos(θ + π + α)
arrow4_y = y2 + length * sin(θ + π + α)
# 矢印を描画
plot!(dx .+ [arrow1_x, x1, arrow2_x], dy .+ [arrow1_y, y1, arrow2_y], color=color, lw=lw, label="")
plot!(dx .+ [arrow3_x, x2, arrow4_x], dy .+ [arrow3_y, y2, arrow4_y], color=color, lw=lw, label="")
# 寸法などを記入
annotate!( (x1 + x2)/2 + dx + deltax, (y1 + y2)/2 + dy + 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, lw=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=lw)
else
plot!(ra .* cosd.(θ) .+ ox, rb .* sind.(θ) .+ oy,
linecolor=color, linestyle=lty, linewidth=lw,
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=lw)
else
plot!(cosine .* x .- sine .* y .+ ox,
sine .* x .+ cosine .* y .+ oy,
linecolor=color, linestyle=lty, linewidth=lw,
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;
##### 正多角形を描く関数
# (ox, oy) を中心，r を半径とする円に内接する正 n 多角形。fcol を指定すると塗りつぶし。
function polygon(ox, oy, r, n; φ=90, color=:black, lty=:solid, lw=0.5, fcolor="")
θ = range(0, 2π, length=n+1) .+ φ / 180 * π
if fcolor == ""
plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
linecolor=color, linestyle=lty, linewidth=lw)
else
plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
linecolor=color, linestyle=lty, linewidth=lw,
seriestype=:shape, fillcolor=fcolor)
end
end;
##### 直線と円の交点座標（2つあることを前提とする）
# 半径が r，中心座標が (x, y) の円と，
# 二点 (x1, y1), (x2, y2) を結ぶ直線が 2 点で交わるとき，
# 交点座標のセット (x3, y3), (x4, y4) を求める関数
function intersection2(r, x, y, x1, y1, x2, y2)
# 直線をベクトル表示
dx = x2 - x1
dy = y2 - y1
# 円の中心を原点に平行移動
fx = x1 - x
fy = y1 - y
# 2次方程式 at² + bt + c = 0 を作る
a = dx^2 + dy^2
b = 2 * (fx * dx + fy * dy)
c = fx^2 + fy^2 - r^2
# 判別式（2点交差前提なので省略してもよい）
D = b^2 - 4a*c
# 解の計算（2点分）
sqrtD = sqrt(D)
t1 = (-b + sqrtD) / (2a)
t2 = (-b - sqrtD) / (2a)
# パラメトリック直線上の点を計算
x3 = x1 + t1 * dx
y3 = y1 + t1 * dy
x4 = x1 + t2 * dx
y4 = y1 + t2 * dy
return (x3, y3, x4, y4)
end;
