/*放物の方程式（抵抗力あり）
  dx/dt=vx
  dz/dt=vz
  dvx/dt=-A/M*vx*sqrt(vx*vx+vz*vz)
  dvz/dt=-A/M*vz*sqrt(vx*vx+vz*vz)-G
  x0=0
  z0=0
  v0x=v0*cos(th)
  v0z=v0*sin(th)
  をルンゲクッタ法で解く
*/
#include <stdio.h>
#include <math.h>

#define TMAX 200.0
#define TH 60
#define V0 35
#define M   1
#define A   0.01
#define G   9.8

double f1(double t,double x,double vx);
double f2(double t,double z,double vz);
double f3(double t,double vx,double vz);
double f4(double t,double vx,double vz);

int main(void)
{
        double x,z,vx,vz,t,dt;
        double k0[4],k1[4],k2[4],k3[4];
        FILE *output;

        dt=0.02;
        x=z=0.0;
        vx=V0*cos(TH*3.1415926/180);
        vz=V0*sin(TH*3.1415926/180);             /* 初期値 */

        output=fopen("deg60.data","w");

        for(t=0.0;t<TMAX;t+=dt) {
                k0[0]=dt*f1(t,x,vx);
                k0[1]=dt*f2(t,z,vz);
                k0[2]=dt*f3(t,vx,vz);
                k0[3]=dt*f4(t,vx,vz);

                k1[0]=dt*f1(t+dt/2.0,x+k0[0]/2.0,vx+k0[2]/2.0);
                k1[1]=dt*f2(t+dt/2.0,z+k0[1]/2.0,vz+k0[3]/2.0);
                k1[2]=dt*f3(t+dt/2.0,vx+k0[2]/2.0,vz+k0[3]/2.0);
                k1[3]=dt*f4(t+dt/2.0,vx+k0[2]/2.0,vz+k0[3]/2.0);

                k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,vx+k1[2]/2.0);
                k2[1]=dt*f2(t+dt/2.0,z+k1[1]/2.0,vz+k1[3]/2.0);
                k2[2]=dt*f3(t+dt/2.0,vx+k1[2]/2.0,vz+k1[3]/2.0);
                k2[3]=dt*f4(t+dt/2.0,vx+k1[2]/2.0,vz+k1[3]/2.0);

                k3[0]=dt*f1(t+dt,x+k2[0],vx+k2[2]);
                k3[1]=dt*f2(t+dt,z+k2[1],vz+k2[3]);
                k3[2]=dt*f3(t+dt,vx+k2[2],vz+k2[3]);
                k3[3]=dt*f4(t+dt,vx+k2[2],vz+k2[3]);

                x=x+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                z=z+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                vx=vx+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                vz=vz+(k0[3]+2.0*k1[3]+2.0*k2[3]+k3[3])/6.0;

                if (z<=0.0) break;

                fprintf(output,"%f %f\n",x,z);

        }

        fclose(output);

        return 0;
}

double f1(double t,double x,double vx)
{
        double r;

        r=vx;

        return r;
}

double f2(double t,double z,double vz)
{
        double r;

        r=vz;

        return r;
}

double f3(double t,double vx,double vz)
{
        double r;

        r=-A/M*vx*sqrt(vx*vx+vz*vz);

        return r;
}

double f4(double t,double vx,double vz)
{
        double r;

        r=-A/M*vz*sqrt(vx*vx+vz*vz)-G;

        return r;
}