/*
 Eliza Dubroff
 ps7#2
 beam.cc
 4/24/97
 */

#include <math.h>
#include "beam.h"

void beam::RK4(double h)
{
  double ku[4];
  double kv[4];
  double u[int((L/h)+1)];
  double v[int((L/h)+1)];
  double x[int((L/h)+1)];
  int i;

  u[0]=0;
  x[0]=0;
  v[0]=0;

for (i=0;i < int(L/h); i++)
  {
    printf("x = %f\tv= %f\n",x[i],(v[i]*12));
    kv[0]=u[i];
    ku[0]=((Q+w*(L-x[i])/2)*(L-x[i])*pow((1+u[i]*u[i]),1.5))/(E*I);
    kv[1]=u[i]+h*ku[0]/2;
    ku[1]=(Q+w*(L-x[i]-h/2)/2)*(L-x[i]-h/2);
    ku[1]=ku[1]*pow((1+(u[i]+h*ku[0]/2)*(u[i]+h*ku[0]/2)),1.5)/(E*I);
    kv[2]=u[i]+h*ku[1]/2;
    ku[2]=(Q+w*(L-x[i]-h/2)/2)*(L-x[i]-h/2);
    ku[2]=ku[2]*pow((1+(u[i]+h*ku[1]/2)*(u[i]+h*ku[1]/2)),1.5)/(E*I);
    kv[3]=u[i]+h*ku[2];
    ku[3]=(Q+w*(L-x[i]-h)/2)*(L-x[i]-h);
    ku[3]=ku[3]*pow((1+(u[i]+h*ku[2])*(u[i]+h*ku[2])),1.5)/(E*I);
    
    x[i+1]=x[i]+h;
    u[i+1]=u[i]+h*(ku[0]+2*ku[1]+2*ku[2]+ku[3])/6;
    v[i+1]=v[i]+h*(kv[0]+2*kv[1]+2*kv[2]+kv[3])/6;   
  }
printf("x = %f\tv= %f\n",x[int((L/h))],(v[int(L/h)]*12));
}