//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// and differentiation (derivative) functions
// Dr. Turhan Coban
// =====================================================
import java.io.*;


class fm1 extends f_xi
{
 double gamma;
	public fm1(double gammai)
    {gamma=gammai;}
 // sürtünmeli eğik atış problemi
 double func(double x[],int x_ref)
 {
   //x[0] is x, x[1] is y
   double g=9.806;
   if(x_ref==0) return x[3];
   else if(x_ref==1) return x[4];
   else if(x_ref==2) return -gamma*Math.sqrt(x[3]*x[3]+x[4]*x[4])*x[3];
   else              return -g - gamma*Math.sqrt(x[3]*x[3]+x[4]*x[4])*x[4];
 }
}

class diferansiyel3F
{
public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {a[j][i]=fi[j-1][i];}
}

return a;
}
public static void main(String args[]) throws IOException
  {
  fm1 b1=new fm1(0.0001);
  double y[];
  y=new double[4];
  double V0=900.0;//  m/s
  double alfa=5.0/180.0*Math.PI; //6 derece
  y[0]=0.0;
  y[1]=0.0;
  y[2]=V0*Math.cos(alfa);
  y[3]=V0*Math.sin(alfa);
  double a[][]=RK6(b1,0.0,200.0,y,10000);
  Plot p1=new Plot(a[1],a[2]);
   double xa[]=new double[2];
  double ya[]=new double[2];
  xa[0]=0;xa[1]=12.0;
  ya[0]=0;ya[1]=0.0;
  p1.addData(xa,ya);
  p1.setXlabel("x");
  p1.setYlabel("y");
  p1.setPlabel("sürtünmeli top atışı "+alfa*180/Math.PI+" derece açı ve "+V0+" m/s mermi çıkış hızı");
  p1.setMinMax(0,12000,0,500);
  p1.plot();
 }
}
