Julia 関数ソースプログラム

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;