HSPTV!掲示板

弾着に関すれうプログラムを書こうとしているのですがそもそもどういう計算式にすれば良いのか、でてこないため困っています

初速はXm/秒
空気抵抗は1秒毎に*0.999
重力加速度はG m/s2

この条件で入力した距離飛ばすには仰角を何度で撃ち出せば良いのか、というものを作りたいです
よろしくお願いいたします

もし、検索のキーワードで困っているのなら
「斜方投射」で検索してみてください。
空気抵抗を考慮した数式ならWikipediaに載ってますね。

斜方投射は真っ先に調べましたが、正直この式をプログラム状の関数や計算式に置き換えることは自分の技術力では難しい...というか不可能でして
上記の条件式自体も理論的にもあまり物理的に現実性があるわけではなくあくまでもそれっぽく動いているように見せるためのもの、で、wikipediaやその他の本格的なサイトに乗っている式の代用日買い物が有ります

Cで書いたやつ

#include <stdio.h>
#include <math.h>

/* 重力加速度 [ m/s/s ] */
#define g 9.8
 /* π */
#define pi 3.141592653589793238462623383279

int main (void)
{
    double t;    /* 経過時間 [ s ]   */
    double v;    /* 初速度   [ m/s ] */
    double x, y; /* 位置     [ m ]   */
    double z;    /* 投射角度 [ rad ] */

    printf("初速度, 投射角度(度), 経過時間");
    scanf("%lf, %lf, %lf", &v, &z, &t);
    z *= pi / 180.0; /* 度 → ラジアン変換 */

    x = v * cos(z) * t;
    y = v * sin(z) * t - g * pow(t, 2.0) / 2.0;

    printf("x = %f\ny = %f\n", x, y);
    return 0;
}

ググって出てきたものですが
一度テストしたいのなら
printf = mes
scanf = 入力
double = .0を変数に
cos sinはHSP標準

まともに数式一撃でキめたいのなら、
プログラムの処理順に合った経過処理回数からなる数列を用意して一般式を解き、方程式に持ち込みましょう。

/*
dx0 = X軸方向初速
dx(n) = dx(n-1)*(1.0-処理時間当たりの空気抵抗率)
dy0 = Y軸方向初速
dy(n) = dy(n-1)*(1.0-処理時間当たりの空気抵抗率)+処理時間当たりの重力加速
　　　※空気抵抗が先か重力加速が先かはプログラムに因る。
x(n) = Σ(i=0->n) dx(i)
y(n) = Σ(i=0->n) dy(i)
　　　※加算末端が 0 or 1 、n or n-1 はプログラムに因る。
y(n),x(n)の一般式が解けたらx(n)=飛距離,y(n)=0の連立式を射出角度について解きましょう。
　　nの候補が複数でるので、0以上の最小の実数値でいいかと思う。
　　整数とか見なかったことにしておきましょう。　　*/

簡単そうに見えて実際に解いてみると値がけっこう扱いづらい形で出てくる上に、
方程式になる手前、代入操作等でうんざりげんなりぐんにょりげっそりします。


と、もはやどんな値なのかよくわからない数式まみれの状態と格闘するのもひとつですが・・・

計算負荷がでかく付きますが二分法の選択もありだと思う。
(練習として角度を変えて何発か打って、本番の角度の決定する感じ)
この場合、解の存在が確実でなければ利用できないので、届かない場合は前もってはじかなければいけません。

最終的にある程度の誤差に収まるような仰角を得られたら計算を打ち切ればよく、わかりやすいんで。
(なお、計算ルーチンは軌道作成の関数を流用できるならしたほうが安全です。)

もちろん計算速度上妥協できない状況ならがんばるしかないです。たぶん。

> この条件で入力した距離飛ばすには仰角を何度で撃ち出せば良いのか、というものを作りたいです

確認です。
発射地点と着弾点は同じ高さですか？
風邪の影響は無視していいですよね？

まともに計算する以外だと、弾道の予想軌道を描画してプレイヤーの視覚に任せるとか、適当に撃ち込んて着弾観測してごにょごにょとかいろいろありそうです。
でもこういうのはやっぱりまともにぶつかりたい。後でやろう。

お正月におこたでぬくぬくしながらスクリプト書いてみました。
そのまえに計算式の説明します。スクリプトは最後です。

------------------------------------------------
変数を次のように定義する
重力加速度：g
t秒後（tフレーム後）の移動距離：x,y
射出時の初速：v0
発射角度：r [rad]

ただし上向きを正とする。
したがって重力加速度gは g<0 である。
着弾点は射出点との高低差0とする。

時刻tのときの距離座標x,yは次式のようになる。

x = vx0 * t
y = 1/2*g*t^2 + vy0 * t

着弾点と射出点は高低差0なので y = 0 のときの t を求める
（なお高低差0以外ならy=高低差として計算すれば良い。）

0 = 1/2*g*t^2 + vy0 * t
= t(1/2*g*t + vy0)

(1/2*g*t + vy0) = 0
1/2*g*t = -vy0
t = -vy0 * 2/g

y = 0となるtは、
t = 0
t = -vy0 * 2/g
t = 0は射出時なので、t = -vy0 * 2/g となる。
このときxは

x = vx0 * t
= vx0 * -vy0 * 2/g

射出角をrとしたとき
vx0 = v0 * cos(r)
vy0 = v0 * sin(r)

したがってxは次のようになる。

x = -v0^2 * cos(r) * sin(r) * 2/g

r が決まれば距離 x を算出することができる。
しかし求めたいものは x が決まっている場合の r なので、この式から r を算出する。

x = -v0^2 * cos(r) * sin(r) * 2/g
cos(r) * sin(r) = x / (-v0^2) * g/2

ここで加法定理
sin(x+y) = sin(x) * cos(y) + cos(x) * sin(y)
y=xの場合
sin(2*x) = 2 * sin(x) * cos(x)
1/2 * sin(2*x) = sin(x) * cos(x)

したがって
cos(r) * sin(r) = x / (-v0^2) * g/2
1/2 * sin(2*r) = x / (-v0^2) * g/2
sin(2*r) = x / (-v0^2) * g
2 * r = asin( x / (-v0^2) * g )
r = 1/2 * asin( x / (-v0^2) * g )
となり角度を求めることができる。r は複数無限に出るがどれを用いても良い。
尚HSP3のasin関数なら解は１個出てくる。

解が出るかどうかの判定は次の式からわかる。
sin(2*r) = x / (-v0^2) * g

-1 <= sin(2*r) <= 1
なので
-1 <= x / (-v0^2) * g <= 1
が満たされれば算出可能で有効射程距離。

-------------------
以上説明終わり。

ということで、スクリプト書いてみたのでスクリプトは長いのでこーどったーに上げてみました。
http://codetter.com/?p=1157

・・・個人的な意見っというかすごいどうでもいい話だと思われますが・・・
高校の現代物理じゃない方の物理を勉強すれば大体の計算式がイメージできると思うのですが・・・

もちろん、HSP用に式を打っていくとしたらまた話は別だとは思いますが・・・


知ってたら申し訳ないですが　物理の運動方程式に関しては基礎で躓くことも多いですが
以外とやって行けば簡単みたいなことも多いと思うので　やってみてはいかがかなーと思います。
年齢とかはわからないので　勝手な推測で描きました。すみません。

空気抵抗を考慮していない回答がちらほら見られるのはなぜでしょうか...。

質問者さんは空気抵抗をどのようにイメージしているのでしょうか？
質問者さんのイメージでは速度ベクトルの反対向きに0.999の加速度が掛かっているのだろうなと想像しましたが、
現実は少し違うと思います。砲弾(丸いと仮定)の速度が十分遅い場合は速度の大きさの定数倍(この定数は砲弾の半径と空気の粘性で決まります。)に比例します。
もちろん向きは速度ベクトルの逆向きで。

この場合は運動方程式が立てられますし、幸運なことに、この方程式においては発射からt秒後の位置を積分で解析的に求められます。
Newton力学を微積分を用いて本格的に議論するのは大学からですが、この運動方程式ならちょっと賢い高校三年生でも解けます。

問題の解決にはなりませんが、私自身の物理の復習も兼ねて、空気抵抗を考慮した場合の発射からt秒後の位置を関数で与えてみました。

http://www.mediafire.com/download/0y2dg2ehndcok5q/%E7%A9%BA%E6%B0%97%E6%8A%B5%E6%8A%97%E3%82%92%E8%80%83%E6%85%AE%E3%81%97%E3%81%9F%E6%96%9C%E6%96%B9%E6%8A%95%E5%B0%84.pdf

解けたつもりですが、間違っている可能性もあります。お気づきの点があればご指摘ください。

この方程式から着弾時間を仰角θの関数で求め、弾体の水平移動距離を所望の値にするようなθを求める(もちろん、そのようなθが存在しない場合もあります。)という操作を行えばよいわけですが、これを解析的に行うのは無理ではないかと思います。
(そういえば高校の頃、これをテーマにして課題学習&発表をやっている友人がいたっけな...。)
条件がそろえば解析解が求まるらしいという話を聞いたことはありますが、私にはわかりません。

結局私としては、堀木さんのアプローチが最も現実的な気がします。

実際のコードは思うに、
A.現在速度を弾の位置に加算
B.現在速度の鉛直成分に重力加速を施す
C.現在速度に空気抵抗による減衰処理を行う(勝手に積算と思ってた)
の3要素をループしていると勝手に想像しています(順番はどうなのか不明ですが)。
GENKIさんもこう考えているとコードを見る限りで分かりました。

なのでしたら、なぜプログラムでの実装方法に則った方法で位置を計算しないのでしょうか。
そうしなければ、推測演算として成立しない気がするんですが。

そして、GENKIさんのコードを見る限り、
あたも当然のように「机上の理屈に合うように実際の動きを補正」しています。
プログラムの在り方ってそういうものなのですか？

私個人としては、
実際の動きのコードがあり、それの予測演算を行うんじゃないの？と首を傾げるばかりです。
(そう思って数列によるアプローチを提示しました)

そう考えると重力加速をy+=(1/2)*g*t^2と置いたり、微分方程式を解いてもその先と結びつかないように思えます。
プログラムでの軌道計算は数列和であって積分じゃないですし、現在速度も微分じゃない。
こういうのって階差微分というんでしたっけ？(どうでもいい)

なお、GENKIさんのコードから補正部分を取っ払うと、重力加速部分はy+=(1/2)*g*t*(t-1)になるはず。
そして数学的に補正が1/2でぴったりなのはさすがです。


なんか批判的な話になってしまい申し訳ないです。
ですが、予測演算の指針としてはっきりしておきたいので。


余談：こちらが提示した数列から解いた場合、代数的に解ける範囲の外のような解が出てきます。
真っ向勝負するとなると絶望的。

個人的にはアクションゲームの敵の挙動に使えるくらいの
処理速度が確保できるかどうかに興味があります。
（たとえば遠く離れた自分めがけて敵が砲弾を撃ってくるときの射出角の計算とか）
60fpsでゲーム全体の処理のループをさせているときに、
発射角の計算で処理落ちすることってないですかね？

>空気抵抗は1秒毎に*0.999
これを完全に見落としていました。(´・ω・`)
またあとでスクリプト考えてみます。

たくさんのご回答、ありがとうございます
nepisatさんと堀木さんの回頭を元に一定の成果が出たら返信を行うと思っていましたが一度まとめて返信させていただきますね
ちなみに現在の状況は堀木さんの2分法のヒントを元にベクトルの分解を使って構築中でした（あんまり進んでいない）
＞堀木さん
計算順序としてはそうなりますね…
t秒後の位置というのはそれなりに求められそうな気もしますが、最初に位置を決めて底に向けて飛ばすというのは単純な逆算では無理そうで
仰るとおり真っ向勝負以外の方法になりそうです

＞GENKIさん
発射点と着弾点の高さは同じものとしての計算を前提に考えています
スクリプトですが非常に参考になります。後は空気抵抗値を入れれば...という感じですがそれがかなり邪魔してるんですよね…

＞alcaneさん
高校では物理は取っていないんですよね…
数学は一応2Bまではやったのでベクトル分解などをの手段をとってみたりと錯誤しているのですが

＞FunnyMakerさん
自分の考えている空気抵抗はその通りです
実際はもっと複雑で、弾道学として学問が成立するレベルなのでそれっぽく見せるための妥協点、としてこのような数値になっています

＞ZAPさん
素人が組むとどの程度の重さになるかはちょっと想像が付きませんね

＞堀木さん

空気抵抗考慮する前のものと同じ要領で計算してみようかと思ったのですが、砲弾の軌道が放物線じゃなくなるのってどうやればいいんだ？と気づいた時点で自力は諦めました。
次にWikipediaの斜方投射の式を解こう！と思ったのですが一見して超めんどくさそうだったので物理の式から直接求めるのはあきらめました。
探せばθを算出した式が見つかるかもしれませんが、掘木さんのレス読んでたらそれもそうだと思って中断です。
正月休みにちょっとやってみるっていうレベルじゃなかった。orz

さて言われて気づきましたがゲームで使うんだろうから厳密解求めたって仕方ないんですよね。
ゲーム内の方が正であり「実際」なので、ゲーム内の結果にあうような解法の方が都合がいい正しいと思います。

> 計算負荷がでかく付きますが二分法の選択もありだと思う。
> (練習として角度を変えて何発か打って、本番の角度の決定する感じ)

ということでこれが一番だと思います。
撃って外れたら補正、を何度か繰り返す。これなら上空で風が吹いてても、障害物に衝突して軌道が変わっても計算できます。（発射前の状態を保存・再現できればObaQでも同じ手が使えるはず。）
何回繰り返しても精度が良くならなければ「届かない」として処理しちゃえばいいと思います。


> 実際の動きのコードがあり、それの予測演算を行うんじゃないの？と首を傾げるばかりです。

そうですよね、普通そうです。
距離から射出角を言い当てるゲーム（弾丸アニメーションを表示しないゲーム）でも作ってるのかなー、変わってるけどまあそういうこともあるかー、とか思ってました。手元で電卓はじいて遊ぶゲームなのかなーとか。
今更ですがそれってどんなゲームだよとわたしも思います。w


> あたも当然のように「机上の理屈に合うように実際の動きを補正」しています。

おまけで付けたアニメーション部分（スクリプトのほぼすべて）と発射角の予測計算は解法が全然違うんで答えもズレます。
なのでおまけのほうを結果に合わせて感で誤魔化そうとしたわけです。
掘木さんが首を傾げるのももっともで、私は物理の計算式の結果を「実際」として答えてしまっています。ゲームのほうが実際と考えるべきでした。


> 砲弾(丸いと仮定)の速度が十分遅い場合は速度の大きさの定数倍(この定数は砲弾の半径と空気の粘性で決まります。)に比例します。

ちょっと気になるので補足。

流体中を移動する物体が流体から受ける抗力は2種類あります。
粘性抵抗と慣性抵抗です。
粘性抵抗は速度に比例して力は大きくなります。
慣性抵抗は速度の二乗に比例して力は大きくなります。
このため速度が遅い場合は粘性抵抗の方の影響が大きく、速い場合は慣性抵抗の影響が大きくなります。
どちらの影響が大きくなるかは計算してみないとわかりません。どちらが効くかわからないなら両方入れちゃえば問題ありません。パラメータ増えるので設定が面倒そうですが…。

しかし空気は粘性が小さいですし砲弾なら速度も速いので、個人的意見としては粘性抵抗は無視できるほど小さくなるんじゃないかと思います。わたしの中のイメージでは粘性抵抗は油中や水中をゆっくり動く時ぐらいしか効かない感じ。あ、計算してないんであくまで個人のイメージです。(^ ^;

HSP開発wikiに参考になりそうな資料がいくつかありました。
http://wiki.hsp.moe/%E7%89%A9%E7%90%86%E3%81%AE%E3%81%8A%E3%81%8A%E3%81%BE%E3%81%8B%E3%81%AA%E8%A9%B1.html
http://wiki.hsp.moe/%E7%89%A9%E7%90%86%E3%81%AE%E3%81%8A%E3%81%8A%E3%81%BE%E3%81%8B%E3%81%AA%E8%A9%B1%EF%BC%8F%E9%A3%9B%E3%81%B6%E3%83%9C%E3%83%BC%E3%83%AB%E3%81%9D%E3%81%AE%EF%BC%93.html
http://wiki.hsp.moe/%E7%89%A9%E7%90%86%E3%81%AE%E3%81%8A%E3%81%8A%E3%81%BE%E3%81%8B%E3%81%AA%E8%A9%B1%EF%BC%8F%E9%A3%9B%E3%81%B6%E3%83%9C%E3%83%BC%E3%83%AB%E3%81%9D%E3%81%AE%EF%BC%94.html
http://wiki.hsp.moe/%E7%89%A9%E7%90%86%E3%81%AE%E3%81%8A%E3%81%8A%E3%81%BE%E3%81%8B%E3%81%AA%E8%A9%B1%EF%BC%8F%E5%B0%91%E3%81%97%E7%9C%9F%E9%9D%A2%E7%9B%AE%E3%81%AA%E8%B3%87%E6%96%99.html
もういくつか消えちゃってますね。一部だけでも消えないようにしてくださってるだけでも個人的には大変ありがたいです。

>ちょっと気になるので補足。〜

GENKIさん、ご指摘ありがとうございます。言われてみればそのとおりですね。
嫌な予感がしたのでちゃんと計算してみたら、砲弾くらい速いと 慣性抵抗 >> 粘性抵抗 でした。

Zaraさんが「それっぽく見せるための妥協点」としているからには多分ゲームを作るつもりなのだろうと私は勝手に思っていますが、
そうであれば私も最初から、細かい理屈をこねてどうのこうの言うよりも空気抵抗が dv/dt = -v*0.999 という設定を尊重すべきだったなと、今更ながら思いました。
回答者が質問者の意向を勝手にねじ曲げて話を進めてはいけませんね...。どうも失礼しました。
ゲームの一部として考えているなら、微小時間毎の速度の積算としての「ゲーム内での弾の位置」に対して正確な予測を出さなければ意味がないですね。

ゲームなら弾の初速も武器に応じて有限種類で良さそうですし、それならば発射地点と着弾地点の高さが有限種類でよいなら限界到達距離は予めシミュレーションで求めて定数として記憶しておけばいいと思います。
これで、弾が届くかどうかを直ぐ判別できますし、数値解の探索に二分法を使うとしても予め求めておいた「最も遠くに飛ぶ仰角θ」と0radの間から始めればいい気がしますね。同じ水平距離飛ばすなら短時間で着弾したほうが嬉しいかなと。(目標の装甲の角度を考慮する等 凝った作りにするならこれじゃダメかも)

数列式いじりまわしてたらそれっぽい値がでる方法に当たったので書き込み
でも残念、これは計算式じゃなくて算出手続きなんだ・・・。

別に二分法じゃなくてもよくね？と適当にやってみたら思った以上に理屈が通ってた。
これはあくまで手続きに過ぎないので、余談として書き込み。

#include "hspmath.as"

#const	SCANNUMINSECOND	60 // 1秒当たりの処理回数

eps = 0.00001				// 計算打ち切り仰角差
g = 2.0 / SCANNUMINSECOND	// 重力加速

p = powf(0.999,1.0/SCANNUMINSECOND)	// 空気抵抗
w = 60.0 / SCANNUMINSECOND	// 初速

h = 0.0				// 目標地点y変位(下方向正)
d = 29.0			// 目標地点x変位(正値でないと動かない)

font msgothic,10

c=0.0	// 初期値(仰角)
n=0.0	// 初期値(到達まで処理回数)

repeat 100	// 無限ループの回避
	// nを割り出す式を水平方向速度から計算
	n=logf((d*p)/(cos(c)*w)-d/(cos(c)*w)+1)/logf(p)
	// 割り出したnで所定の高さに合う角度を計算
	sc = -( (h*(1.0-p)-g*n)/(1.0-powf(p,n)) + g/(1.0-p) )/w
	prevc = c
	c = asin(sc)
	mes "仰角"+c+"/到達まで"+n
	// 角度が非数の時点でエラー
	if ( isnan(c) ) : break
	// 角度の差が一定以下になったら打ち切り
	if absf(prevc - c ) < eps : break
loop

// 実際に飛ばして確認する部分
pos 180,0
dy = -w*sin(c)
dx = w*cos(c)		// 初速の定義
repeat 120,1
	h -= dy			// 位置のシフト
	d -= dx
	dy *= p			// 空気抵抗
	dx *= p
	dy += g			// 重力加速
	tmpstr = strf("%2d",cnt)+strf("|目標まで (%9.5f,%9.5f)",d,h)
	mes tmpstr+strf("//速度(%9.5f,%9.5f)",dx,dy)
loop

//	　手続きの理屈上、真横(仰角0)で打ちだした際に目標地点の上を通る場合に
//	正しく計算できないことが確定している。
//	　この制約に当たる状況は、1処理時間の初速よりも目標地点が近い場合である。
//	(つまりは近すぎる場合)
//	この場合は何かしらの例外処理で補完する必要がある。
//	　また、オマケ程度で目標地点の高さも指定できるようにしたが、
//	下方向に打ちおろす場合に特に制約に引っかかる事象が多く、
//	役に立たないのが現実。
//	
//	また、double型特有の演算誤差の影響で正しい値が得られないこともあるし、
//	　精度を上げ過ぎると失敗するのはこの手の計算の定番といえるか。

遠くに届くかどうかはほぼ確実に検出できるようになったものの、余裕で届きすぎる場合がダメとは
これが灯台下暗しというやつか(違