sagantaf

メモレベルの技術記事を書くブログ。

円錐曲線(楕円、双曲線、放物線)のグラフをgo言語で描く

背景

go言語と宇宙工学の勉強のために、円錐曲線のグラフを描いてみました。

円錐曲線の式で一気に、放物線も双曲線も楕円も描けます。

$$ r = \frac{l}{1 + e cos\theta} $$

離心率eと半直弦lを変えることで、放物線、双曲線、楕円、円になります。

main関数

まずは、イメージを湧かせるためにmain()関数から記載します。

func main(){
    // データ取得
    thetaList := GetThetaData()

    // パラメタ指定
    e := 0.5 // 離心率
    l := 1.5 // 半直弦
    graphTitle = "parabola polar"
    savePath = "/tmp/parabola-polar.png"

    // 描画実行
    DrawConic(e, l, thetaList, graphTitle, savePath)
}

流れとしては、thetaの数値リストを取得し、パラメタを指定した後に、描画用の関数を実行しています。

thetaデータの生成

ここでは単純に0° ~ 360°の時のthetaを10°ごとに計算し、リストにしているだけです。

func GetThetaData() (thetaData []float64){
    for angle := 0.0; angle <= 360; {
        theta := angle * math.Pi / 180
        thetaData = append(thetaData, theta)
        angle = angle + 10
    }
    return
}

DrawConic()関数

円錐曲線を計算するクロージャCalcConicPolar()(後述)を設定後、データ1つ1つに対して適用していきます。

また、極座標を直交座標に変換しています。

最後に自作関数plotDrawLiner()(後述)にて描画とpng保存を実行します。

func DrawConic(e, l float64, thetaList []float64, graphTitle, savePath string){
    // クロージャを設定
    c := CalcConicPolar(e, l)

    // radius算出
    var radiusList []float64
    for _, v := range thetaList{
        r := c(v)
        radiusList = append(radiusList, r)
    }

    // 座標変換
    var xys [][]float64
    for i, theta := range thetaList{
        radius := radiusList[i]
        xy := make([]float64, 2)
        x, y := Polar2cartesian(radius, theta)
        xy[0] = math.Round(x*100) / 100
        xy[1] = math.Round(y*100) / 100
        xys = append(xys, xy)
    }

    // 描画&保存
    PlotDrawLiner(graphTitle, savePath, xys)
}

クロージャ CalcConicPolar()関数

極座標による円錐曲線を計算するためのクロージャ(極座標θ, r)です。

離心率: e, 半直弦: l をクロージャの引数として受け取るようになっています。

このクロージャを呼び出すごとにthetaを渡すことで、radiusが返されます。

func CalcConicPolar(e, l float64) func(theta float64) float64 {
    return func(theta float64) float64 {
        return l / (1 + e * math.Cos(theta)) // 冒頭の円錐曲線の式より
    }
}

座標変換(極座標→直交座標)Polar2cartesian()関数

グラフにプロットできるようにするために、極座標r, θから直交座標x, yへ変換する関数をつくります。

func Polar2cartesian(r, theta float64)(x, y float64){
    x = r * math.Cos(theta)
    y = r * math.Sin(theta)
    return
}

グラフの描画と保存 PlotDrawLiner()関数

グラフの描画には、

GitHub - gonum/plot: A repository for plotting and visualizing data

を利用しました。

func PlotDrawLiner(graphTitle, savePath string, data [][]float64){
    // 図の作成
    p, err := plot.New()
    if err != nil{
        log.Fatalln("error: ", err)
    }

    // グラフの装飾の設定
    p.Title.Text = graphTitle
    p.X.Label.Text = "X"
    p.Y.Label.Text = "Y"
    p.Add(plotter.NewGrid()) // 補助線

    // 描画するデータの取得
    pts := make(plotter.XYs, len(data))
    var xDataList []float64
    var yDataList []float64

    for i, pair := range data {
        pts[i].X = pair[0]
        xDataList = append(xDataList, pair[0])
        pts[i].Y = pair[1]
        yDataList = append(yDataList, pair[1])
    }

    // 折れ線グラフの描画
    err = plotutil.AddLinePoints(p, pts)
    if err != nil {
        log.Fatalln("failed to draw liner: ", err)
    }

    // 座標範囲の設定
    xMax, xMin, err := GetMaxMin(xDataList)
    if err != nil {
        log.Fatalln("failed to get max and min from xDataList: ", err)
    }
    yMax, yMin, err := GetMaxMin(yDataList)
    if err != nil {
        log.Fatalln("failed to get max and min from yDataList: ", err)
    }
    p.X.Min = xMin - 1
    p.X.Max = xMax + 1
    p.Y.Min = yMin - 1
    p.Y.Max = yMax + 1

    // 保存
    if err := p.Save(6*vg.Inch, 6*vg.Inch, savePath); err != nil {
        log.Fatalln("Failed to save plot:", err)
    }
}

最大値と最小値を取得 GetMaxMin()関数

与えられたリストの中で最大値と最小値を返す関数です。

標準ライブラリには、二つの数値を比較する関数しか存在しなかったため、自作しました。

ちなみにその標準ライブラリを中で利用しています。(math.Max, math.Min)

func GetMaxMin(inputList []float64)(max, min float64, err error){
    if len(inputList) == 0 {
        return 0, 0, errors.New("input list is empty")
    }
    max = inputList[0]
    min = inputList[0]
    for _, p := range inputList{
        max = math.Max(max, p)
        min = math.Min(min, p)
    }
    err = nil
    return
}

コード全体

最後にここまでのコード全体を記載します。

package main

import (
    "gonum.org/v1/plot"
    "gonum.org/v1/plot/plotter"
    "gonum.org/v1/plot/plotutil"
    "gonum.org/v1/plot/vg"
    "log"
    "math"
)

func GetThetaData() (thetaData []float64){
    for angle := 0.0; angle <= 360; {
        theta := angle * math.Pi / 180
        thetaData = append(thetaData, theta)
        angle = angle + 10
    }
    return
}

func CalcConicPolar(e, l float64) func(theta float64) float64 {
    return func(theta float64) float64 {
        return l / (1 + e * math.Cos(theta)) // 冒頭の円錐曲線の式より
    }
}

func PlotDrawLiner(graphTitle, savePath string, data [][]float64){
    // 図の作成
    p, err := plot.New()
    if err != nil{
        log.Fatalln("error: ", err)
    }

    // グラフの装飾の設定
    p.Title.Text = graphTitle
    p.X.Label.Text = "X"
    p.Y.Label.Text = "Y"
    p.Add(plotter.NewGrid()) // 補助線

    // 描画するデータの取得
    pts := make(plotter.XYs, len(data))
    var xDataList []float64
    var yDataList []float64

    for i, pair := range data {
        pts[i].X = pair[0]
        xDataList = append(xDataList, pair[0])
        pts[i].Y = pair[1]
        yDataList = append(yDataList, pair[1])
    }

    // 折れ線グラフの描画
    err = plotutil.AddLinePoints(p, pts)
    if err != nil {
        log.Fatalln("failed to draw liner: ", err)
    }

    // 座標範囲の設定
    xMax, xMin, err := GetMaxMin(xDataList) // common.goの関数
    if err != nil {
        log.Fatalln("failed to get max and min from xDataList: ", err)
    }
    yMax, yMin, err := GetMaxMin(yDataList) // common.goの関数
    if err != nil {
        log.Fatalln("failed to get max and min from yDataList: ", err)
    }
    p.X.Min = xMin - 1
    p.X.Max = xMax + 1
    p.Y.Min = yMin - 1
    p.Y.Max = yMax + 1

    // 保存
    if err := p.Save(6*vg.Inch, 6*vg.Inch, savePath); err != nil {
        log.Fatalln("Failed to save plot:", err)
    }
}


func DrawConic(e, l float64, thetaList []float64, graphTitle, savePath string){
    c := CalcConicPolar(e, l)

    // radius算出
    var radiusList []float64
    for _, v := range thetaList{
        r := c(v)
        radiusList = append(radiusList, r)
    }

    // 座標変換
    var xys [][]float64
    for i, theta := range thetaList{
        radius := radiusList[i]
        xy := make([]float64, 2)
        x, y := Polar2cartesian(radius, theta)
        xy[0] = math.Round(x*100) / 100
        xy[1] = math.Round(y*100) / 100
        xys = append(xys, xy)
    }

    // 描画&保存
    PlotDrawLiner(graphTitle, savePath, xys)
}


func main(){
    // データ取得
    thetaList := GetThetaData()

    // 楕円
    e := 0.5 // 離心率
    l := 1.5 // 半直弦
    graphTitle = "parabola polar"
    savePath = "/tmp/parabola-polar.png"

    // 描画実行
    DrawConic(e, l, thetaList, graphTitle, savePath)
}

結果

これを実行することで、以下のようなpngファイルを取得できます。

f:id:sagantaf:20200824004134p:plain
楕円

また、e=1.1など、離心率eを1以上にすることで、双曲線が描けます。

f:id:sagantaf:20200824004104p:plain
双曲線