// Simulation des Bremsvorganges eines Verkehrsflugzeugs
// http://walter.bislins.ch/doku/aviatik

function CTurbofan() {
  // data sheet
  this.F_max    = 117900;  // N
  this.m0_max   = 397.3;   // m/s
  this.mu_max   = 6;
  // parameters from other sources
  this.F_idle   = 3500;    // N
  this.theta    = 60.0 / 180.0 * Math.PI;
  this.v9_max   = 343;     // m/s = a0
  this.N1_max   = 1.0;
  this.N1_idle  = 0.215;
  this.N1_maxrev= 0.8;
  this.N2_max   = 0.972;
  this.N2_idle  = 0.607;
  // input parameters
  this.N1       = this.N1_idle;
  this.v0       = 0.0;
  // simulated values as functions of N1 and v0
  this.F        = 0.0;    // N
  this.F_rev    = 0.0;    // N
  // simulated internal values as functions of N1 and v0
  this.N2       = 0.0;
  this.F1       = 0.0;    // N
  this.F2       = 0.0;    // N
  this.m1       = 0.0;    // kg/s
  this.m2       = 0.0;    // ks/s
  this.v9       = 0.0;    // m/s
  this.v19      = 0.0;    // m/s
  this.mu       = 0.0;
  this.kv       = 0.0;
  // internal constant parameters, computed in Init()
  this.k        = 0.0;
  this.r        = 0.0;
  this.k1       = 0.0;
  this.r1       = 0.0;
  this.a        = 0.0;
  this.b        = 0.0;
  this.km1      = 0.0;
  this.km2      = 0.0;
  this.m1_max   = 0.0;
  this.m1_idle  = 0.0;
  this.m2_max   = 0.0;
  this.kv_max   = 0.0;
  this.kv_idle  = 0.0;
  this.F1_max   = 0.0;
  this.F1_idle  = 0.0;
  this.v9_idle  = 0.0;
}

CTurbofan.prototype.Init = function(){
  this.m1_max = this.m0_max / (this.mu_max + 1);
  this.m2_max = this.m1_max * this.mu_max;
  this.a = (this.N2_max - this.N2_idle) / (this.N1_max - this.N1_idle);
  this.b = this.N2_idle - this.a * this.N1_idle;
  this.k = this.SolveForK( this.N1_max, 1.0, this.N1_idle, this.F_idle/this.F_max );
  this.r = this.N1_max / Math.log( (1.0 / this.k) + 1.0 );
  this.km1 = this.m0_max / ((this.mu_max + 1.0) * (this.a + this.b));
  this.km2 = this.m0_max * this.mu_max / (this.mu_max + 1.0);
  this.kv_max = ((this.F_max * (1.0 + this.mu_max)) / (this.m0_max * this.mu_max * this.v9_max)) - (1.0 / this.mu_max);
  this.kv_idle = this.KVofN1( this.N1_idle );
  this.F1_max = this.m1_max * this.v9_max;
  this.m1_idle = this.M1ofN1( this.N1_idle );
  this.v9_idle = this.FofN1(this.N1_idle) / ( this.M1ofN1(this.N1_idle) + this.M2ofN1(this.N1_idle) * this.kv_idle );
  this.F1_idle = this.m1_idle * this.v9_idle;
  this.k1 = this.SolveForK( this.N2_max, 1.0, this.N2_idle, this.F1_idle/this.F1_max );
  this.r1 = this.N2_max / Math.log( (1.0 / this.k1) + 1.0 );
}

CTurbofan.prototype.Compute = function( n1, v ) {
  if (xDef(n1)) { this.N1 = n1; }
  if (xDef(v)) { this.v0 = v; }
  this.N2 = this.N2ofN1(this.N1);
  this.F1 = this.F1ofN1(this.N1);
  this.F2 = this.F2ofN1(this.N1);
  this.m1 = this.M1ofN1(this.N1);
  this.m2 = this.M2ofN1(this.N1);
  this.v9 = this.F1 / this.m1;
  this.v19 = (this.m2 < 0.01) ? 0.0 : (this.F2 / this.m2);
  this.mu = this.m2 / this.m1;
  this.kv = this.v19 / this.v9;
  this.F  = this.m1 * (this.v9 - this.v0) + this.m2 * (this.v19 - this.v0);
  this.F_rev = this.m1 * (this.v9 - this.v0) - this.m2 * (Math.cos(this.theta) * this.v19 + this.v0);
}

CTurbofan.prototype.FofN1 = function( n ) {
  // require k, r, is computed
  return this.k * (Math.exp(n/this.r) - 1.0) * this.F_max;
}

CTurbofan.prototype.F1ofN1 = function( n ) {
  // require k1, r1, F1_max is computed
  return this.k1 * (Math.exp(this.N2ofN1(n)/this.r1) - 1.0) * this.F1_max;
}

CTurbofan.prototype.F2ofN1 = function( n ) {
  // require k, r, k1, r1, F1_max is computed
  return this.FofN1(n) - this.F1ofN1(n);
}

CTurbofan.prototype.N2ofN1 = function( n ) {
  return this.a * n + this.b;
}

CTurbofan.prototype.M1ofN1 = function( n ) {
  // require km1 is computed!
  return this.km1 * this.N2ofN1(n);
}

CTurbofan.prototype.M2ofN1 = function( n ) {
  // require km2 is computed!
  return this.km2 * n;
}

CTurbofan.prototype.KVofN1 = function( n ) {
  // require kv_max is computed!
  return (this.kv_max - 1.0) * n + 1.0;
}

CTurbofan.prototype.SolveForK = function( N1, F1, N2, F2 ) {
  var myFunction = function(x) {
    var t1 = N1 * Math.log( (F2 / x) + 1 );
    var t2 = N2 * Math.log( (F1 / x) + 1 );
    return t1 - t2;
  };
  var myErrFunction = function( aMessage ) {
    this.Status = aMessage;
  }
  this.Status = '';
  var k = SolveWithNewton( myFunction, 0.0000001, 0.0000001, myErrFunction );
  if (this.Status == '') {
    this.Status = 'ok';
  } else {
    alert( this.Status );
  }
  return k;
}

function CNewtonSolver() {
  this.Reset();
}

CNewtonSolver.prototype.Reset = function() {
  this.F = function(X) { return X; }
  this.Guess = 0.0;
  this.Precision = 0.0001;
  this.Eps = this.Precision / 2.0;
  this.Result = 0.0;
  this.MaxLoops = 1000;
  this.Status = 'not solved';
  this.NLoops = 0;
}

CNewtonSolver.prototype.Solve = function( aFunction, aGuess, aPrecision ) {
  var dx, X, Y, Xnew, df;
  if (typeof(aFunction)  !='undefined') { this.F         = aFunction;  }
  if (typeof(aGuess)     !='undefined') { this.Guess     = aGuess;     }
  if (typeof(aPrecision) !='undefined') { this.Precision = aPrecision; }
  this.NLoops = 0;
  this.Status = '';
  this.Eps = this.Precision / 2.0;
  Xnew = this.Guess;
  for (dx = this.Precision*2.0; (dx > this.Precision) && (this.NLoops < this.MaxLoops); dx = Math.abs( Xnew - X )) {
    X = Xnew;
    try {
      Y = this.F(X);
      df = (this.F(X + this.Eps) - Y) / this.Eps;
      Xnew = X - (Y / df);
    } catch(err) {
      this.Status = err.message;
      return this.Result;
    }
    this.NLoops++;
  }
  this.Result = Xnew;
  if (this.NLoops >= this.MaxLoops) {
    this.Status = 'max loops exceedet';
  } else {
    this.Status = 'solved';
  }
  return this.Result;
}

function SolveWithNewton( aFunction, aGuess, aPrecision, aErrFunction ) {
  // function aFunction( X )
  // function aErrFunction( aMessage )
  var solver = new CNewtonSolver();
  solver.Solve( aFunction, aGuess, aPrecision );
  if (solver.Status != 'solved') {
    if (aErrFunction) { aErrFunction( solver.Status ); }
  }
  return solver.Result;
}