/* Copyright (C) 2001-2008 Claudio Girardi
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


/* This is a collection of various formulas for the inductance of 
 * loops and coils
 */

/* The following formulas are from
 *   R. Lundin, 
 *   "A Handbook Formula for the Inductance of a Single-Layer Circular Coil", 
 *   Proc. IEEE, vol. 73, no. 9, pp. 1428-1429, Sep. 1985
 */

function f1(x) {
  return (1 + x * (0.383901 + 0.017108 * x)) / (1 + 0.258952 * x);
}


function f2(x) {
  return (x * (0.093842 + x * (0.002029 - x * 0.000801)));
}

/* IndCalc(a, b, n) - returns the inductance of a multiturn single layer circular coil
 * a: coil radius [m] 
 * b: coil length [m] 
 * n: number of turns 
 */
function IndCalc(a, b, n) {
  var u0, shape_ratio, l;

  u0 = 0.4 * Math.PI * 1e-06;
  shape_ratio = 2 * a / b;

  if (shape_ratio <= 1) {
    l = u0 * n * n * Math.PI * a * a * (f1(shape_ratio * shape_ratio) - (4 /  (3 * Math.PI)) * shape_ratio) / b;
  } else {
    l = u0 * n * n * a * ((Math.log(4 * shape_ratio) - 0.5) * f1(1 / (shape_ratio * shape_ratio)) + f2(1 / (shape_ratio * shape_ratio)));
  }
  return l;
}


/* The following formulas are from
 *   F. E. Terman,
 *   "Radio Engineers' Handbook", London,
 *   McGraw-Hill, 1st edition, 1950
 */

/* MutIndConc(a1, a2, b1, b2, n1, n2) - return the mutual inductance between two coaxial concentrical multiturn single layer coils
 * a1: larger radius coil radius [m] 
 * b1: larger radius coil length [m] 
 * n1: larger radius coil number of turns 
 * a2: smaller radius coil radius [m] 
 * b2: smaller radius coil length [m] 
 * n2: smaller radius coil number of turns 
 */
function MutIndConc(a1, a2, b1, b2, n1, n2) {
  /* The formula was converted to metric units */
  /* A term was added, in order to match the coaxial coils */
  /* inductance formula (see below) */

  if (b2 > b1) {
    /* exchange b2 and b1 */
    var b_tmp = b1;
    b1 = b2;
    b2 = b_tmp;
  }

  var g = Math.sqrt(a1 * a1 + b1 * b1 / 4);

  var M = 1.972441e-06 * a2 * a2 * n1 * n2 * (1 + (a1 * a1 * a2 * a2 / (8 * g * g * g * g)) * (3 - b2 * b2 / (a2 * a2)) +

 (a1 * a1 * a1 * a1 * a2 * a2 * a2 * a2) * (3 - b1 * b1 / (a1 * a1)) * (2.5 - 2.5 * b2 * b2 / (a2 * a2) + b2 * b2 * b2 * b2 / (4 * a2 * a2 * a2 * a2)) / (32 * g * g * g * g * g * g * g * g)

) / g;

  return M;
}

/* MutIndCoax(a1, a2, b1, b2, n1, n2, D) - return the mutual inductance between two coaxial non-concentrical multiturn single layer coils
 * a1: larger radius coil radius [m] 
 * b1: larger radius coil length [m] 
 * n1: larger radius coil number of turns 
 * a2: smaller radius coil radius [m] 
 * b2: smaller radius coil length [m] 
 * n2: smaller radius coil number of turns 
 * D: axial distance between centers of coils [m]
 */
function MutIndCoax(a1, a2, b1, b2, n1, n2, D) {
  /* The formula was converted to metric units */
  /* A couple of typos in the original formula (?) */
  /* were corrected */

  if (b2 > b1) {
    /* exchange b2 and b1 */
    var b_tmp = b1;
    b1 = b2;
    b2 = b_tmp;
  }

  var x1 = D - b1 / 2;
  var x2 = D + b1 / 2;
  var r1 = Math.sqrt(x1 * x1 + a1 * a1);
  var r2 = Math.sqrt(x2 * x2 + a1 * a1);

  var K1 = (2 / (a1 * a1)) * (x2 / r2 - x1 / r1);
  var k1 = b2;

  var K3 = -0.5 * (x1 / (r1 * r1 * r1 * r1 * r1) - x2 / (r2 * r2 * r2 * r2 * r2));
  var k3 = 0.5 * a2 * a2 * b2 * (3 - b2 * b2 / (a2 * a2));

  var K5 = -0.125 * a1 * a1 * ((x1 / (r1 * r1 * r1 * r1 * r1 * r1 * r1 * r1 * r1)) * (3 - 4 * x1 * x1  / (a1 * a1)) - (x2 / (r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2 * r2)) * (3 - 4 * x2 * x2  / (a1 * a1)));
  var k5 = 0.5 * a2 * a2 * a2 * a2 * b2 * (2.5 - 2.5 * b2 * b2 / (a2 * a2) + b2 * b2 * b2 * b2 / (4 * a2 * a2 * a2 * a2));

  var M = 0.9862205e-06 * a1 * a1 * a2 * a2 * n1 * n2 * (K1 * k1 + K3 * k3 + K5 * k5) / (b1 * b2);

  return M;
}


/* The following formula is from
 *   F. Langford-Smith (ed.),
 *   "Radiotron Designer's Handbook", 4th edition
 *   Wireless Press, 1952.
 */

function CoilQ(a, b, f) {
  var Q = Math.sqrt(f) / (0.069 / a + 0.054 / b);

  return Q;
}


/* RectLoopIndRoundWire(a, b, d) - returns the inductance of a single turn rectangular loop of round wire
 * a: rectangle side length [m]
 * b: rectangle other side length [m]
 * d: wire diameter [m]
 */
function RectLoopIndRoundWire(a, b, d) {
  /* the formula was converted to metric units */
  var g, l;

  g = Math.sqrt(a * a + b * b); 
  l = 0.4 * ((a + b) * Math.log(4.0 * a * b / d) - a * Math.log(a + g) - b * Math.log(b + g)) + 0.4 * (2.0 * g + d - 2.0 * (a + b));

  return l;
}

/* RectLoopIndRectWire(a, b, s1, s2) - returns the inductance of a single turn rectangular loop of rectangular wire
 * a: rectangle side length [m]
 * b: rectangle other side length [m]
 * s1: wire thickness [m]
 * s2: wire width [m]
 */
function RectLoopIndRectWire(a, b, s1, s2) {
  /* the formula was converted to metric units */
  var g, l;
  
  g = Math.sqrt(a * a + b * b);
  l = 0.4 * ((a + b) * Math.log(2.0 * a * b / (s1 + s2)) - a * Math.log(a + g) - b * Math.log(b + g)) + 0.4 * (2.0 * g + 0.447 * (s1 + s2) - (a + b) / 2.0);

  return l;
}


/* RegLoopInd(p, d, idx) - returns the inductance of a single turn loop of regular shape
 * p: loop perimeter [m]
 * d: loop wire diameter [m]
 * idx: loop type index (see below)
 */
function RegLoopInd(p, d, idx) {
  var theta, l;

  switch (idx) {
    case 0: 
      theta = 2.451; /* circle */
    break;
    case 1:  
      theta = 2.561; /* regular octagon */
    break;
    case 2:
      theta = 2.636; /* regular hexagon */
    break;
    case 3:
      theta = 2.712; /* regular pentagon */
    break;
    case 4:
      theta = 2.853; /* square */
    break;
    case 5:
      theta = 3.197; /* equilateral triangle */
    break;
    case 6:
      theta = 3.332; /* isosceles right-angled triangle */
    break;
    default:
      theta = 100; /* error */
   } 
    
   l = 2.0e-07 * p * (Math.log(4.0 * p / d) - theta);

   return l;
}

function toPrec(x) {
  return x.toFixed(3 - Math.log(x)/Math.LN10);
}