ｋ 離心率ω0 最下点角速度
ｌ 半径
ｇ 重力加速度
ｖ０ 最下点速度
θ 最下点位置から ｔ 秒後の角度
ローラー コースター
ローラーコースターの速度と、回転の計算
ここでは、単純に理論的な事についてのみ計算します。
移動体は、慣性モーメントを持たない、点質量です。
侵入速度が遅い場合は、単振り子の周期の計算と同じになりますので、そちらも参照してください。
ｋ 離心率
ω0 最下点角速度
ｌ 半径
ｇ 重力加速度
ｖ０ 最下点速度
θ 最下点位置から ｔ 秒後の角度
離心率ｋの値によって、往復するか回転するかが変わります。
最下点の速度が、コースターの直径の高さから落下した時の速度と同じ場合、
ｋ＝１となり、最上点で留まります。
ｋ＜１の場合は単振り子と同じになり、往復運動をします。
ｋ＞１の時は、最上点を超えて、回転運動となります。
ｋ＜１の単振り子の場合の時間Ｔは、往復してもとへ戻る時間で、ｋ＞１の場合は、一回転の時間です。
ｋ＝１の時は最上点に達するには、無限大の時間を要することを表しています。
左図は、離心率ｋの値が１を超えず往復している場合です。
赤丸は、最下点通過後の経過時間位置を示しています。
離心率ｋが１になるように値を設定するのは難しいので、チェックを入れることによって、離心率ｋを１にして計算することが出来ます。
運動の様子は、アニメーションで確認できますが、時間は、計算値より長くなります。
プログラム
楕円積分が第一二種完全楕円積分となっていますが、第一種の値しか使用されません。
//----------------------------------------------------- // ヤコビの楕円関数 Roller Coaster // http://fuchino.ddo.jp/yatsugatake/ellipticx.pdf // http://www.wannyan.net/scidog/ellipse/ellipse.pdf //----------------------------------------------------- unit Main; interface uses Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.ExtCtrls, VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, VCLTee.TeeProcs, VCLTee.Chart; type TForm1 = class(TForm) LabeledEdit1: TLabeledEdit; LabeledEdit2: TLabeledEdit; Memo1: TMemo; TimeBtn: TButton; LabeledEdit3: TLabeledEdit; Image1: TImage; LabeledEdit4: TLabeledEdit; StartAnimeBtn: TButton; EndAnimeBtn: TButton; Timer1: TTimer; Label1: TLabel; CheckBox1: TCheckBox; procedure TimeBtnClick(Sender: TObject); procedure FormCreate(Sender: TObject); procedure LabeledEdit1Change(Sender: TObject); procedure LabeledEdit2Change(Sender: TObject); procedure LabeledEdit4Change(Sender: TObject); procedure StartAnimeBtnClick(Sender: TObject); procedure EndAnimeBtnClick(Sender: TObject); procedure Timer1Timer(Sender: TObject); private { Private 宣言 } function jacobi_sn(u, k: double): double; procedure Complete_Elliptic_Integrals(k: double; var Fk, Ek: double); procedure CanvasClear; procedure DrawCircle(x, y, r: double; w: integer; c: Tcolor; var xi, yi: integer); procedure DrawLine(xs, ys, xe, ye: double; w: integer; c : Tcolor); function calc(L, V, t, k: double; var LT: double; var Ch: integer): double; procedure textOut(x, y: double; s: string); procedure TimeCircle(x, y, r: integer); public { Public 宣言 } end; var Form1: TForm1; implementation {$R *.dfm} uses Math; var DBL_EPSILON : double; // 桁落ち判定値 InOutXY : array of array of integer; // 直線区間座標配列 CircleXY : array of array of integer; // 円ループ座標配列 MF : byte; // 0 停止 1 start 2 Loop 3 end 4 タイマー停止 CON : integer; // 配列カウント iolength : integer; // 直線区間配列長さ looplength : integer; // ループ区間配列長さ xb, yb : integer; // 消去用image座標 ERASEF : boolean; // 消去フラグ true 消去 StopF : boolean; // ループ脱出フラグ true 脱出 Q, Qs : double; // Q 角度 Qs 指定時間角度 xlb : double; // 直線部分x位置 kF : boolean; // k > 1 true k < 1 false // 第一二種完全楕円積分 // K 離心率 // Fk 第一種 // Ek 第二種 procedure TForm1.Complete_Elliptic_Integrals(k: double; var Fk, Ek: double); const KN = 30; var m : double; // the parameter of the elliptic function m = k^2 a : double; // arithmetic mean g : double; // geometric mean a_old : double; // previous arithmetic mean g_old : double; // previous geometric mean two_n : integer; // power of 2 sum : double; half : double; // 0.5 N : integer; // loop数 begin half := 1 / 2; if k = 0 then begin Fk := PI / 2; Ek := Fk; exit; end; m := k * k; if m = 1 then begin Fk := infinity; Ek := 1; exit; end; a := 1; g := sqrt(1 - m); two_n := 1; sum := (2 - m); N := 0; while N < KN do begin // 収束しなかった時KN回で終了 g_old := g; a_old := a; a := (g_old + a_old) * half; g := sqrt(g_old * a_old); two_n := two_n + two_n; sum := sum - two_n * (a * a - g * g); if abs(a_old - g_old) <= a_old * DBL_EPSILON then break; inc(N); end; Fk := PI / 2 / a; Ek := PI / 4 / a * sum; end; // Jacobi関数 Jacobi sn(u, k) function TForm1.Jacobi_sn(u, k: double): double; const KN = 30; var a: array of double; g: array of double; c: array of double; two_n : integer; phi : double; half : double; i, j : integer; begin half := 1 / 2; if k = 0 then begin result := sin(u); exit; end; if k = 1 then begin result := sin(2 * arctan(exp(u)) - pi / 2); exit; end; setlength(a, KN + 1); setlength(g, KN + 1); setlength(c, KN + 1); a[0] := 1; g[0] := sqrt(1 - k * k); c[0] := k; two_n := 1; for i := 0 to KN do begin // 収束しなかった時KNで終了 if abs(a[i] - g[i]) < a[i] * DBL_EPSILON then break; two_n := two_n + two_n; a[i + 1] := half * (a[i] + g[i]); g[i + 1] := sqrt(a[i] * g[i]); c[i + 1] := half * (a[i] - g[i]); end; phi := two_n * a[i] * u; for j := i downto 1 do phi := half * (phi + arcsin(c[j] * sin(phi) / a[j])); result := sin(phi); end; // イメージ全消去 procedure TForm1.CanvasClear; begin Image1.Canvas.Brush.Style := bsSolid; Image1.Canvas.Brush.Color := clWhite; // Image1.Canvas.Brush.Color := clYellow; Image1.Canvas.FillRect(Rect(0,0,Image1.Width,Image1.Height)); end; // 円作図 // x,y 原点座標 // r 半径 // w 線の幅 // c 線の色 // xi, yi image座標 procedure TForm1.DrawCircle(x, y, r: double; w: integer; c : Tcolor; var xi, yi: integer); var HC, VC : integer; ri : integer; begin HC := Image1.Width div 2; VC := Image1.Height div 2; xi := round(x) + HC; yi := VC - round(y); ri := round(r); image1.Canvas.pen.Style := psSolid; image1.Canvas.Pen.Color := c; image1.Canvas.Pen.Width := w; image1.Canvas.Ellipse(xi - ri, yi - ri, xi + ri + 1, yi + ri + 1); end; // 直線作図 // xs, ys 始点原点座標 // xe, ye 終点原点座標 // w 線の幅 // c 線の色 procedure TForm1.DrawLine(xs, ys, xe, ye: double; w: integer; c : Tcolor); var HC, VC : integer; xsi, ysi : integer; xei, yei : integer; begin HC := Image1.Width div 2; VC := Image1.Height div 2; xsi := round(xs) + HC; ysi := VC - round(ys); xei := round(xe) + HC; yei := VC - round(ye); image1.Canvas.pen.Style := psSolid; image1.Canvas.Pen.Color := c; image1.Canvas.Pen.Width := w; Image1.Canvas.MoveTo(xsi, ysi); Image1.Canvas.LineTo(xei, yei); end; // 文字出力 // x,y 原点座標 // s 文字 procedure TForm1.textOut(x, y: double; s: string); var HC, VC : integer; xi, yi : integer; begin Image1.Canvas.Font.Size := 12; HC := Image1.Width div 2; VC := Image1.Height div 2; xi := round(x) + HC - 5; yi := VC - round(y) - 5; image1.Canvas.TextOut(xi, yi, s); end; // 時間と角度計算 // L 半径 // V 初速度 // t 経過時間 // k 離心率 // LT 周期時間 // ch 回転数 又は 周期数 function TForm1.calc(L, V, t, k: double; var LT: double; var Ch: integer): double; const g = 9.80665; // 重力加速度 m/s^2 var Fk, Ek : double; // 第一二種楕円積分値 u, Q : double; begin kF := False; // 回転フラグリセット result := 0; if k = 1 then begin // 頂点で停止 LT := infinity; u := sqrt(g / L) * t; Q := k * Jacobi_sn(u, k); result := arcsin(Q) * 2 / pi * 180; end; if k < 1 then begin // 往復移動 Complete_Elliptic_Integrals(k, Fk, Ek); // 第一二種楕円積分 LT := 4 * sqrt(L / g) * Fk; ch := trunc(t / LT); u := sqrt(g / L) * t; Q := k * Jacobi_sn(u, k); result := arcsin(Q) * 2 / pi * 180; end; if K > 1 then begin // 回転移動 Complete_Elliptic_Integrals(1 / k, Fk, Ek); // 第一二種楕円積分 LT := 2 / k * sqrt(L / g) * FK; ch := trunc(t / LT); t := t - LT * ch; u := sqrt(g / L) * k * t; Q := Jacobi_sn(u, 1 / k); Q := arcsin(Q) * 2 / pi * 180; result := Q; if t > LT / 2 then result := 360 - Q; kF := True; // 回転フラグセット end; textOut(-2, -125, '0°'); // 角度表示 textOut(112, 5, '90°'); if K >= 1 then begin textOut(-10, 130, '180°'); textOut(-130, 5, '270°'); end else begin textOut(-15, 130, '±180°'); textOut(-130, 5, '-180°'); end; end; // 各種計算 procedure TForm1.TimeBtnClick(Sender: TObject); const g = 9.80665; // 重力加速度 m/s^2 var k : double; // 離心率 ch, i : integer; L : double; // 半径 V : double; // 初速度 t : double; // 経過時間 LT : double; // 周期時間 st, dt : double; // 時間計算用 dt = Δt x, y : double; // 原点0,0の座標 Np : integer; // 分割数 xi, yi : integer; // Image座標 begin MF := 0; ERASEF := false; startAnimeBtn.Enabled := false; endAnimeBtn.Enabled := false; timer1.Enabled := false; val(Labelededit1.Text, L, ch); if ch <> 0 then begin application.MessageBox('コースターの半径に間違いがあります。','注意',0); exit; end; if L <= 0 then begin application.MessageBox('コースターの半径はゼロ以上にして下さい。','注意',0); exit; end; val(Labelededit2.Text, V, ch); if ch <> 0 then begin application.MessageBox('初速度の値に間違いがあります。','注意',0); exit; end; val(Labelededit3.Text, t, ch); if ch <> 0 then begin application.MessageBox('経過時間の値に間違いがあります。','注意',0); exit; end; val(Labelededit4.Text, Np, ch); if ch <> 0 then begin application.MessageBox('作図分割数の値に間違いがあります。','注意',0); exit; end; if Np < 10 then begin application.MessageBox('作図分割数のが小さすぎますす。','注意',0); exit; end; if Np mod 2 <> 0 then begin application.MessageBox('作図分割数は偶数にしてください。','注意',0); exit; end; image1.Canvas.Pen.Mode := pmCopy; CanvasClear; memo1.Clear; Ch := 0; Qs := 0; LT := 0; // 離心率 if checkbox1.Checked = true then begin // K=1指定時 V := 2 * sqrt(g * L); // 最下点速度計算 k := 1; // 計算だと誤差で1にならないので1セット memo1.Lines.Add('最下点速度 V = ' + FloattostrF(V, fffixed, 14, 10)); end else k := V / (2 * sqrt(g * L)); memo1.Lines.Add('離心率 k = ' + FloattostrF(k, fffixed, 10, 6)); // 離心率と経過時間の角度 回転数 周期時間 又は 周期数 計算 if k = 1 then begin // k=1 頂点で停止 Qs := calc(L, V, t, k, LT, ch); memo1.Lines.Add('離心率 k = 1'); end; if k < 1 then begin // k<1 往復移動 Qs := calc(L, V, t, k, LT, ch); memo1.Lines.Add('離心率 k < 1'); end; if K > 1 then begin // k>1 回転移動 Qs := calc(L, V, t, k, LT, ch); memo1.Lines.Add('離心率 k > 1'); end; memo1.Lines.Add('一周期時間 T = ' + FloattostrF(LT, fffixed, 10, 6)); memo1.Lines.Add(''); memo1.Lines.Add('経過時間の位置 θ = ' + FloattostrF(Qs, fffixed, 10, 6) + '°'); if K > 1 then memo1.Lines.Add('回転数 = ' + IntTostr(ch)); if K < 1 then memo1.Lines.Add('往復数 = ' + IntTostr(ch)); DrawCircle(0, 0, 150, 2, clGreen, xi, yi); // ループ作図 Image1.Canvas.Brush.Style := bsClear; // 頂点で停止する場合 if LT = infinity then begin x := sin(Qs / 180 * pi) * 150; y := - cos(Qs / 180 * pi) * 150; DrawCircle(x, y, 12, 2, clRed, xi, yi); V := V - V * 1E-8; // 作図の為速度を僅かに小さく k := V / (2 * sqrt(g * L)); // 作図用離心率設定 Qs := calc(L, V, t, k, LT, ch); // LT作図用用周期計算 dt := LT / Np; // 作図時間間隔 for i := 0 to Np div 4 do begin // 一周期の4分の一作図 st := dt * i; Q := calc(L, V, st, k, LT, ch); // 移動体作図と角度取得 x := sin(Q / 180 * pi) * 150; y := - cos(Q / 180 * pi) * 150; DrawCircle(x, y, 10, 1, clBlack, xi, yi); // 移動体作図 end; i := 0; repeat // 直線部分の移動体作図 inc(i); x := 150 / L * V * dt * i; DrawCircle(-x,-150, 10, 1, clBlack, xi, yi); // 移動体作図 until x > 150; DrawLine(-X, -150, X, -150, 2, clBlue); // 直線作図 exit; end; // ループ部分移動イメージ作図 と 座標配列作成 dt := LT / Np; // 作図時間間隔 // 直線部分 座標配列長さ設定 i := 0; repeat inc(i); x := 150 / L * V * dt * i; until x > 150; if i < 2 then begin application.MessageBox('作図分割数のが小さすぎますす。','注意',0); exit; end; setlength(InOutXY, i, 2); // 直線部分配列確保 iolength := i - 1; // 配列添え字は0から xlb := x; timer1.Interval := round(dt * 1000); // timer時間設定msec setlength(CircleXY, Np + 1, 2); // ループ配列長さ確保 Looplength := NP; for i := 0 to Np do begin st := dt * i; // 経過時間 Q := calc(L, V, st, k, LT, ch); // 移動体作図と角度取得 x := sin(Q / 180 * pi) * 150; y := - cos(Q / 180 * pi) * 150; DrawCircle(x, y, 10, 1, clBlack, xi, yi); // 移動体作図 CircleXY[i, 0] := xi; // 移動体座標 x CircleXY[i, 1] := yi; // 移動体座標 y end; // 経過時間の位置の赤〇作図 x := sin(Qs / 180 * pi) * 150; y := - cos(Qs / 180 * pi) * 150; DrawCircle(x, y, 12, 2, clRed, xi, yi); memo1.Lines.Add('赤丸 経過時間の位置'); DrawLine(-xlb, -150, xlb, -150, 2, clBlue); // 直線作図 // 直線部分移動イメージ作図 と 座標配列 i := 0; repeat inc(i); x := 150 / L * V * dt * i; DrawCircle(-x,-150, 10, 1, clBlack, xi, yi); // 移動体作図 DrawCircle(x,-150, 10, 1, clBlack, xi, yi); // 移動体作図 InOutXY[i - 1, 0] := xi; // 移動体座標 x InOutXY[i - 1, 1] := yi; // 移動体座標 y until x > 150; startAnimeBtn.Enabled := true; CON := 0; end; // 移動体作図 procedure TForm1.TimeCircle(x, y, r: integer); begin image1.Canvas.Ellipse(x - r, y - r, x + r + 1, y + r + 1); end; // タイマーによる移動体作図 アニメーション // フラグMFにより順次作図 procedure TForm1.Timer1Timer(Sender: TObject); begin if ERASEF then begin // 前の円消去 ERASEF := False; TimeCircle(xb, yb, 10); end; if MF = 4 then begin // 終了処理 startAnimeBtn.Enabled := true; MF := 0; Timer1.Enabled := False; end; if (MF = 3) and (CON <= iolength) then begin // 終了直線部分 xb := InOutXY[CON, 0]; yb := InOutXY[CON, 1]; TimeCircle(xb , yb, 10); inc(CON); ERASEF := True; if CON > iolength then MF := 4; // 終了設定 end; if (MF = 2) and (CON <= Looplength) then begin // ループ部分 xb := CircleXY[CON, 0]; yb := CircleXY[CON, 1]; TimeCircle(xb, yb, 10); inc(CON); ERASEF := True; if STOPF then begin // ルーフ脱出処理 if CON > Looplength then begin MF := 3; CON := 0; end; end else if CON = Looplength then CON := 0; end; if (MF = 1) and (CON <= iolength) then begin // スタート直線部分 xb := Image1.Width - InOutXY[iolength - CON, 0]; yb := InOutXY[iolength - CON, 1]; TimeCircle(xb, yb, 10); inc(CON); ERASEF := True; if CON > iolength then begin // ループイン設定 MF := 2; CON := 0; end; end; end; // text変更入力によるアニメーションボタン設定 procedure TForm1.LabeledEdit1Change(Sender: TObject); begin startAnimeBtn.Enabled := false; endAnimeBtn.Enabled := false; end; procedure TForm1.LabeledEdit2Change(Sender: TObject); begin startAnimeBtn.Enabled := false; endAnimeBtn.Enabled := false; end; procedure TForm1.LabeledEdit4Change(Sender: TObject); begin startAnimeBtn.Enabled := false; endAnimeBtn.Enabled := false; end; // アニメーションスタート procedure TForm1.StartAnimeBtnClick(Sender: TObject); var xi, yi : integer; x, y : double; begin startAnimeBtn.Enabled := false; image1.Canvas.Pen.Mode := pmCopy; CanvasClear; Image1.Canvas.Brush.Style := bsClear; DrawCircle(0, 0, 150, 1, clGreen, xi, yi); // ループ円作図 DrawLine(-xlb, -150, xlb, -150, 1, clBlue); // 直線部分描画 x := sin(Qs / 180 * pi) * 150; y := - cos(Qs / 180 * pi) * 150; DrawCircle(x, y, 12, 1, clRed, xi, yi); // 指定経過時間位置円 textOut(-2, -125, '0°'); // 角度表示 textOut(112, 5, '90°'); if KF then begin textOut(-10, 130, '180°'); textOut(-130, 5, '270°') end else begin textOut(-20, 130, '±180°'); textOut(-130, 5, '-90°'); end; textOut(-180, 210, 'アニメーションの時間は計算時間より長くなります。'); CON := 0; image1.Canvas.Pen.Color := clBlack; // 移動体の色黒 image1.Canvas.Pen.Width := 2; image1.Canvas.Pen.Mode := pmNotXor; endAnimeBtn.Enabled := true; MF := 1; // アニメーションフラグ start 1 STOPF := False; // ループ脱出フラグリセット Timer1.Enabled := true; // アニメーションタイマースタート end; // アニメーション停止 procedure TForm1.EndAnimeBtnClick(Sender: TObject); begin endAnimeBtn.Enabled := false; STOPF := true; // ループ停止フラグセット end; // 初期設定 procedure TForm1.FormCreate(Sender: TObject); var bitmap : Tbitmap; y, z : double; i : integer; begin startAnimeBtn.Enabled := false; endAnimeBtn.Enabled := false; timer1.Enabled := false; // 桁落ち判定値設定 i := 0; z := 1; DBL_EPSILON := 1; repeat DBL_EPSILON := DBL_EPSILON / 2; y := z + DBL_EPSILON; inc(i); until (y = 1) or (i > 70); // 桁落ちしたら終了 DBL_EPSILON := DBL_EPSILON * 2; // 桁落ち判定値 // Imageにビットマップ設定 無くても可 bitmap := Tbitmap.Create; bitmap.Width := Image1.Width; bitmap.Height := Image1.Height; Image1.Picture.Bitmap := bitmap; bitmap.Free; TimeBtnClick(nil); end; end.
Roller_coaster.zip