//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// Runge-Kutta 6 method with multible differential equations
// Dr. Turhan Coban
// =====================================================
// dy1/dx=f(x,y1,y2,..,yn)
// dy2/dx=f(x,y1,y2,..,yn)
// ......
// dyn/dx=f(x,y1,y2,..,yn)
// Pendulum problem

import java.io.*;

class fm1 extends f_xi
{
// 
double ca     = 1.55;
double cv     = 519;
double R      = 1.05;
double r      = 0.068;
double Vstr   = 67.9;
double alpha0 = 93;
double alphas = 93;
double alphap = 93;
double alphaH = 0.84;
double beta0  = 7;
double betas  = 7;
double betap  = 7;
double betaH  = 1.17;
double gammaH = 0;
double P0      = 93;
double Paval   = P0;
double Pvval   = (1 / (1 + p.R/p.r)) * P0;
double Hval    = (1 / (p.R * p.Vstr)) * (1 / (1 + p.r/p.R)) * P0;
double tau     = 4;

double func(double x[],int x_ref)
{
if (x[0] <= 600) R = 1.05;
else R = 0.21 * exp(600-x[0]) + 0.84; 
double y;   
if(x_ref==0) 
{y=- (1 / (ca * R)) * x[2] + (1/(ca * R)) * x[2]+ (1/ca) * Vstr * x[3];}
else if(x_ref==1) 
{y=(1 / (cv * R)) * Paoft- ( 1 / (cv * R) + 1 / (cv * r) ) * x[2];}
else if(x_ref==2) 
{
double Ts = 1 / ( 1 + (x[1] / alphas)^betas );
double Tp = 1 / ( 1 + (alphap / Paoft)^betap );
y = (alphaH * Ts) / (1 + gammaH * Tp) - betaH * Tp;
}
 }
}

class NA101A
{

public static void main(String args[]) throws IOException
{
  fm1 b1=new fm1();
  double y[];
  y=new double[3];
  y[0]=Math.PI/16.0; //11.25 degree swing angle
  y[1]=0.0;
  double a[][]=NA93.RK6(b1,0.0,1000.0,y,1000000);
  int n=a[0].length;
  double teta[]=new double[n];
  for(int i=0;i<n;i++) {teta[i]=y[0]*Math.cos(Math.sqrt(b1.g/b1.l)*a[0][i]);}
  Plot p1=new Plot(a[0],a[1]);
  p1.addData(a[0],teta);
  p1.setXlabel("time");
  p1.setYlabel("pendulum angle");
  p1.setPlabel("Pendulum model");
  p1.plot();
 }
}
