◆◆◆ダブルスピン総合スレッド　Part2◆◆◆

Visual C++のコンソールアプリケーションで、空気抵抗とマグナス力を様々に変えて
プログラム作った。当然数値を出力するだけ。0.5秒を50個に区切り、
その間のｘ座標（マウンドからの距離）、y座標（高さ。とりあえず2mにしておいた）
を求めた。発射角は水平で、実際のピッチングと違うが。

「ベースボールの物理学」と言う本を見たら、野球のボール程度の規模だと
空気抵抗は慣性抵抗で速度の2乗に比例するのが普通だが、その本のグラフ
を見たら速度の1次関数か、比例くらいが妥当かと思ったから勝手に比例式に
してしまった（ｽﾏｿ
他にも、マグナス力の係数も粘性抵抗係数も何回か数値入れて適当に
改ざんしてしまった（0.5m先で15cmお辞儀して、）

初速は40m/sをx方向にのみとした。いろんなページや本でプログラム見たら
速度にcos掛けて解いているところが多いが。
まずプログラムのコードを載せる。変なところあったら推敲など頼む。
それから順々に、ジャイロとバックスピンの角度を変えて数値を載せていく。
後でそれ見て、誰かが数値を配列にしてJavaアプレットで
アニメーションにしてくれるかも知れない

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

void runge_kutta(double (*)(double,double,double,double,double),
double (*)(double,double,double,double,double),
double (*)(double,double,double,double,double),
double (*)(double,double,double,double,double),
double,double,double,double,double,double,
double * ,double *,double *,double *);
double func1(double,double,double,double,double);
double func2(double,double,double,double,double);
double func3(double,double,double,double,double);
double func4(double,double,double,double,double);


void main()
{
int i,n;
double h,t,t0,tn,x,xx,y,yy,vx,vxx,vy,vyy;
char file_name [30] ="in.dat";
FILE *fp;

fp=fopen(file_name,"w");
if(fp==NULL){
printf("ファイルをオープンできません\n");
exit(1);
}

printf("===2階の常微分方程式の解法===\n");
printf("---諸元入力---\n");

printf(" t0= ");
scanf("%lf",&t0);
printf(" x0= ");
scanf("%lf",&x);
printf(" y0= ");
scanf("%lf",&y);
printf(" vx0= ");
scanf("%lf",&vx);
printf(" vy0= ");
scanf("%lf",&vy);
printf(" tn= ");
scanf("%lf",&tn);
printf(" n= ");
scanf("%d",&n);

h=(tn-t0)/n;

printf("---解析結果---\n");
printf("t x y");
for(i=0;i<=n;i++)
{
t=t0+i*h;
fprintf(fp," %12.8lf %12.8lf %12.8lf %12.8lf \n",t,x,vx,y);
runge_kutta(func1,func2,func3,func4,t,x,y,vx,vy,h,&xx,&yy,&vxx,&vyy);
x=xx;
y=yy;
vx=vxx;
vy=vyy;
}
fclose(fp);
}

void runge_kutta(double (*Func1)(double,double,double,double,double),
double (*Func2)(double,double,double,double,double),
double (*Func3)(double,double,double,double,double),
double (*Func4)(double,double,double,double,double),
double t,double x,double y,double vx,double vy,double h,
double *xx,double *yy,double *vxx,double *vyy)

{
double k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;

k11=(*Func1)(t,x,y,vx,vy)*h;
k21=(*Func2)(t,x,y,vx,vy)*h;
k31=(*Func3)(t,x,y,vx,vy)*h;
k41=(*Func4)(t,x,y,vx,vy)*h;

k12=(*Func1)(t+h/2.0,x+k11/2.0,y+k21/2.0,vx+k31/2.0,vy+k41/2.0)*h;
k22=(*Func2)(t+h/2.0,x+k11/2.0,y+k21/2.0,vx+k31/2.0,vy+k41/2.0)*h;
k32=(*Func3)(t+h/2.0,x+k11/2.0,y+k21/2.0,vx+k31/2.0,vy+k41/2.0)*h;
k42=(*Func4)(t+h/2.0,x+k11/2.0,y+k21/2.0,vx+k31/2.0,vy+k41/2.0)*h;


k13=(*Func1)(t+h/2.0,x+k12/2.0,y+k22/2.0,vx+k32/2.0,vy+k42/2.0)*h;
k23=(*Func2)(t+h/2.0,x+k12/2.0,y+k22/2.0,vx+k32/2.0,vy+k42/2.0)*h;
k33=(*Func3)(t+h/2.0,x+k12/2.0,y+k22/2.0,vx+k32/2.0,vy+k42/2.0)*h;
k43=(*Func4)(t+h/2.0,x+k12/2.0,y+k22/2.0,vx+k32/2.0,vy+k42/2.0)*h;


k14=(*Func1)(t+h,x+k13,y+k23,vx+k33,vy+k43)*h;
k24=(*Func2)(t+h,x+k13,y+k23,vx+k33,vy+k43)*h;
k34=(*Func3)(t+h,x+k13,y+k23,vx+k33,vy+k43)*h;
k44=(*Func4)(t+h,x+k13,y+k23,vx+k33,vy+k43)*h;