///////////////////////////////////////////////////////////////////////////
//  GeoTC.                                                               //
///////////////////////////////////////////////////////////////////////////
//
//  A library of functions for TC reactance calculations, working from the
//  geometry of the resonator.
//
//  Author:  Paul Nicholson, tssp0807@abelian.org
//
//  The master source for this library is maintained at
//     http://abelian.org/tssp/geotc/
//
//
//  Entry points are prefixed XX_ for clarity, and internal routines and
//  variables are prefixed xx_ to preserve name space.
//
//  Tunable parameters.
//

var xx_VERSION = "2.96";                              // Version of library
var xx_DEBUG = false;                             //  Set true for messages

var GMRES_MAXITER = 100;                  // Max number of gmres iterations 
var GMRES_TOL = 0.001;                          // gmres tolerance target

// 
// Global variables and fixed constants.
//

var EPSILON = 8.854188e-12;                        // Permittivity of space
var MU = (4 * 3.14159265 * 1.0E-7);                // Permeability of space

var xx_errwin = null;                          // Window for error messages
var xx_open_errwin = null;    // External function to open the error window

///////////////////////////////////////////////////////////////////////////
//  Inductance Stuff                                                     //
///////////////////////////////////////////////////////////////////////////
//
//  Tables from F.W. Grover, 'Inductance Calculations', 1946;
//  Reprinted by Dover, 1962.

var grover13 = new Array(    
   0.008297,  
   0.007810, 0.007371, 0.006974, 0.006611, 0.006278,
   0.005970, 0.005685, 0.005420, 0.005173,
   0.004941, 
   0.004723, 0.004518, 0.004325, 0.004142, 0.003969,
   0.003805, 0.003649, 0.003500, 0.003359,
   0.003224, 
   0.003095, 0.002971, 0.002853, 0.002740, 0.0026317,
   0.0025276, 0.0024276, 0.0023315, 0.0022391,
   0.0021502, 
   0.0020646, 0.0019821, 0.0019026, 0.0018259, 0.0017519,
   0.0016805, 0.0016116, 0.0015451, 0.0014808,
   0.0014186, 
   0.0013585, 0.0013004, 0.0012443, 0.0011900, 0.0011374,
   0.0010865, 0.0010373, 0.0009897, 0.0009436,
   0.0008990, 
   0.0008558, 0.0008141, 0.0007736, 0.0007345, 0.0006966,
   0.0006600, 0.0006246, 0.0005903, 0.0005571,
   0.0005251, 
   0.0004941, 0.0004642, 0.0004353, 0.0004074, 0.0003805,
   0.0003545, 0.0003295, 0.0003054, 0.0002823,
   0.00025998, 
   0.00023859, 0.00021806, 0.00019840, 0.00017959, 0.00016162,
   0.00014450, 0.00012821, 0.00011276, 0.00009815,
   0.00008438   );

var grover14 = new Array( 
   0.079093, 
   0.077647, 0.076200, 0.074753, 0.073306, 0.071860,
   0.070413, 0.068966, 0.067520, 0.066073,
   0.064626, 
   0.063180, 0.061733, 0.060287, 0.058840, 0.057394,
   0.055947, 0.054500, 0.053055, 0.051609,
   0.050163, 
   0.048717, 0.047272, 0.045827, 0.044382, 0.042938,
   0.041494, 0.040051, 0.038608, 0.037167,
   0.035727, 
   0.034288, 0.032851, 0.031416, 0.029984, 0.028554,
   0.027128, 0.025707, 0.024291, 0.022881,
   0.021478, 
   0.020084, 0.018700, 0.017329, 0.015972, 0.014632,
   0.013311, 0.012013, 0.010742, 0.009502,
   0.008297  );

var grover15 = new Array(
   .39227 - 9, 
   .54228 - 9, .69229 - 9, .85230 - 9, .99232 - 9, .14234 - 8,
   .29237 - 8, .44240 - 8, .59244 - 8, .74250 - 8,
   .89257 - 8, 
   .04265 - 7, .19276 - 7, .34289 - 7, .49306 - 7, .64327 - 7,
   .79354 - 7, .94388 - 7, .09430 - 6, .22484 - 6,
   .39551 - 6, 
   .54636 - 6, .69744 - 6, .84879 - 6, .00051 - 5, .15268 - 5,
   .30542 - 5, .45891 - 5, .61334 - 5, .76899 - 5,
   .92622 - 5  );

//
//  Compute the mutual inductance between two circular filaments, perpendicular
//  distance d apart, with radii r1 and r2. Our variable K stands for Grover's
//  k' squared.
//

function xx_filament( r1, r2, d, w)
{
   var K = ((r1 - r2)*(r1 - r2) + d*d) / ((r1+r2)*(r1+r2) + d*d);
   var fk, g1, g2;

   if( K < 1e-8)
      return MU * r1 * (Math.log(8*r1/w) - 2) + (MU * r1/4.0);

   if( K < 0.1)
   {
      var i = Math.floor(Math.log(K)*Math.LOG10E * 10);

      if( i <= -60) fk = 0.014468 * ( Math.log(1/K)*Math.LOG10E -0.53307);
      else
      {
         g1 = grover14[ i+60];   g2 = grover14[ i+60-1];
         fk = ((Math.log(K)*Math.LOG10E-(i/10.0))/-0.1) * (g2-g1) + g1;
      }
   }
   else
   if( K >= 0.9)
   {
      var i = Math.floor(Math.log(1-K)*Math.LOG10E * 10);
      if( i <= -40)
         fk = Math.pow( 10.0, .39224 - 3.0 + 1.5 * Math.log(1-K)*Math.LOG10E);
      else
      {
         g1 = grover15[ i+40];   g2 = grover15[ i+40-1];
         fk = Math.pow( 10.0, ((Math.log(1-K)*Math.LOG10E-(i/10.0))/-0.1) * 
                             (g2-g1) + g1);
      }
   }
   else
   {
      var i = Math.floor(K * 100);
      g1 = grover13[ i-10];   g2 = grover13[ i-10+1];
      fk = (K-(i/100.0))/0.01 * (g2-g1) + g1;
   }

   //  The factor 100.0 converts from Grover's cgs to SI, and 1e-6 converts
   //  from uH to H.

   return 1e-6 * 100.0 * Math.sqrt(r1 * r2) * fk;
}

//
//  Compute the self inductance of coil c.
//

function compute_Ldc_thin( c, N, rt)
{
   var t_height = c.height1 + (c.height2 - c.height1) * c.tap/c.turns;
   var t_radius = c.radius1 + (c.radius2 - c.radius1) * c.tap/c.turns;

   var dH = (t_height - c.height1)/N;
   var dR = (t_radius - c.radius1)/N;
   var L = 0;

   for( var j=0; j<N; j++)
   {
      if( this.progress != null) this.progress( Math.floor( 100*j/N));
      for( var k=0; k<N; k++)
      {
         var d = (j - k) * dH;
         var r1 = (j + 0.5) * dR + c.radius1 + rt;
         var r2 = (k + 0.5) * dR + c.radius1 + rt;

         L += xx_filament( r1, r2, d, c.wirad);
      }
   }

   return L * c.tap * c.tap/(N*N);
}

//
//  Compute the self inductance of a ribbon coil c.  The ribbon is
//  decomposed into c.NF parallel filaments spaced evenly across
//  the width.   We solve for the branch currents, then find the total
//  current for 1 volt across the coil.
//

function compute_Ldc_ribbon( c, N)
{
   var rstep = c.rwidth/c.NF;
   var vstep = c.vwidth/c.NF;

   var LM = new Array( N+1);
   for( var j=0; j<c.NF; j++) LM[j+1] = new Array( N+1);

   for( var j=0; j<c.NF; j++)
      for( var k=j; k<c.NF; k++)
         if( j == k)
         {
            var rt = (j+0.5) * rstep;
            var t = compute_Ldc_thin( c, N, rt);
            LM[j+1][j+1] = t;
         }
         else
         {
            var rt1 = (j+0.5) * rstep;
            var rt2 = (k+0.5) * rstep;
            var vt1 = (j+0.5) * vstep;
            var vt2 = (k+0.5) * vstep;
            var t = compute_Mdc_thin( c, c, N, rt1, vt1, rt2, vt2);
            LM[j+1][k+1] = t;
            LM[k+1][j+1] = t;
         }

   var a = new Array( c.NF+1);  // Current through each branch a[1...NF]
   var v = new Array( c.NF+1);  // Voltage across each branch v[1..NF]
   for( j=1; j<=c.NF; j++) v[j] = 1.0;  // 1 volt across the parallel network

   xx_gmres( LM, a, v, c.NF);

   // a[1..NF] now contains the current through each branch 
   var at = 0;
   for( j=1; j<=c.NF; j++) at += a[j];

   var mj = 0.0;
   for( j=1; j<=c.NF; j++)
   {
      a[j] /= at;
      mj += a[j] * j;
   }
   mj /= c.NF;

   c.roffset = c.rwidth * mj;
   c.voffset = c.vwidth * mj;
   c.idist = a;          // Save current distribution

   return 1.0/at;
}

function XX_compute_Ldc( c)
{
   var N = this.ND * this.detail_val;

   if( c.coil == this.COIL_RIBBON)
      return compute_Ldc_ribbon( c, N);
   else
      return compute_Ldc_thin( c, N, 0.0);
}

//
//  Compute the mutual inductance between thin filament coils c1 and c2.
//

function compute_Mdc_thin( c1, c2, N, rt1, vt1, rt2, vt2)
{
   var M = 0;
   var t1_height = c1.height1 + (c1.height2 - c1.height1) * c1.tap/c1.turns;
   var t1_radius = c1.radius1 + (c1.radius2 - c1.radius1) * c1.tap/c1.turns;
   var t2_height = c2.height1 + (c2.height2 - c2.height1) * c2.tap/c2.turns;
   var t2_radius = c2.radius1 + (c2.radius2 - c2.radius1) * c2.tap/c2.turns;

   var c1_dH = (t1_height - c1.height1)/N;
   var c1_dR = (t1_radius - c1.radius1)/N;
   var c2_dH = (t2_height - c2.height1)/N;
   var c2_dR = (t2_radius - c2.radius1)/N;

   for( var j=0; j<N; j++)
   {
      if( this.progress != null) this.progress( Math.floor( 100*j/N));
      for( var k=0; k<N; k++)
      {
         var d = c1.height1 + vt1 + (j + 0.5) * c1_dH - 
                 c2.height1 - vt2 - (k + 0.5) * c2_dH;

         var r1 = (j + 0.5) * c1_dR + c1.radius1 + rt1;
         var r2 = (k + 0.5) * c2_dR + c2.radius1 + rt2;

         M += xx_filament( r1, r2, d, (c1.wirad+c2.wirad)/2);
      }
   }

   return M * c1.tap * c2.tap/(N*N);
}

function XX_compute_Mdc( c1, c2)
{
   var N = this.ND * this.detail_val;
   return compute_Mdc_thin( c1, c2, N, 
                            c1.roffset, c1.voffset,
                            c2.roffset, c2.voffset);
}

///////////////////////////////////////////////////////////////////////////
//  Functions dealing with system components                             //
///////////////////////////////////////////////////////////////////////////
//
//  Functions to specify the resonator geometry.  Each function breaks a
//  resonator component into a number of charge rings.
//
//  These functions can be called several times if you have more than
//  one of each component.
//

//  Install a component p into the system.  We just keep an array system
//  with all the bits in.  Each component is given a unique serial so that
//  we can find it and remove it later, if necessary.

function xx_create_component( geotc, type, nrings)
{
   p = new Object();
   p.type = type;
   p.nrings = nrings;
   p.serial = geotc.serial++;
   geotc.system[ geotc.count++] = p;

   return p;
}

//  Install a coil.  This specifies the winding only, for the purposes of
//  inductance calculations.  To compute capacitances, instead you must
//  call XX_install_primary or XX_install_secondary.

function XX_install_coil( radius1, height1, radius2, height2, turns, wirad)
{
   var c = xx_create_component( this, this.TYPE_COIL, 0);
   c.radius1 = radius1;
   c.radius2 = radius2;
   c.height1 = height1;
   c.height2 = height2;
   c.turns = turns;
   c.tap = turns;
   c.wirad = wirad;
   c.coil = this.COIL_THIN;
   c.NF = 1;
   c.setup = null;
   c.roffset = 0.0;
   c.voffset = 0.0;
  
   return c;
}

//  Install a coil wound from flat ribbon.  rwidth and vwidth specify the
//  radial and vertical projections of the ribbon's width.  rthick is the
//  thickness of the ribbon.

function XX_install_ribbon( radius1, height1, radius2, height2, turns, tap,
                            rwidth, vwidth, rthick)
{
   var c = this.install_coil( radius1, height1, radius2, height2, 
                              turns, rthick);

   c.tap = tap;
   c.rwidth = rwidth;
   c.vwidth = vwidth;
   c.coil = this.COIL_RIBBON;
   c.NF = 5;

   return c;
}

//  Decompose a groundplane into charge rings.

function xx_setup_groundplane( geotc)
{
   var N1 = Math.floor( xx_set_detail( geotc, this)/2);
   var N2 = N1;

   var rm = (this.sec_radius != null) ? this.sec_radius*1.5 : 0;
   if( rm > 0 && rm < this.radius)
   {
      var rm = this.sec_radius * 1.5;
      var hr = rm/N1;

      //  Create N1 rings extending from the center to 1.5 times the
      //  secondary radius.

      for( var k=0; k < N1; k++) 
         new xx_create_ring( geotc, (k+0.5)*hr, 0, hr, 0, 0);

      //  From 1.5 times the secondary radius outwards, create N2 rings.

      hr = (this.radius - rm)/N2;
      for( var k=0; k<N2; k++) 
         new xx_create_ring( geotc, rm + (k+0.5)*hr, 0, hr, 0, 0);
   }
   else
   {
      var N = N1 + N2;
      hr = this.radius/N;
      for( var k=0; k<N; k++) new xx_create_ring( geotc, (k+0.5)*hr, 0, hr, 0, 0);
   }
}

//  Specify a groundplane.  This is slightly different to a grounded disc,
//  because we use a different ring decomposition in order to concentrate
//  rings in the region beneath the secondary (if there is one. If there's
//  no secondary, groundplane degenerates to a grounded disc).

function XX_install_groundplane( radius)
{
   var c = xx_create_component( this, this.TYPE_GND,
                                      this.NG_GND);
   c.radius = radius;
   c.setup = xx_setup_groundplane;

   return c;
}

//  Construct a graded pattern of ring widths to concentrate smaller rings
//  near the ends.

function xx_choose_ring_pattern( N)
{
   var S = new Array( N);
   for( var n = 0; n < N; n++) S[n] = 1/N;

   for( var k = 1; k<=2; k++)
   {
      var e = 0;
      for( var n = 0; n<k; n++)
      {
         S[n] /= 2;
         S[N-1-n] /= 2;
         e += S[n] + S[N-1-n];
      }

      for( var n = k; n < N-k; n++) S[n] += e/(N-2*k);
   }

   return S;
}

//  Decompose a secondary into charge rings.  Uses the same voltage
//  profiles as E-Tesla611.C, which may only be valid for cylinder coils.

function xx_setup_secondary( geotc)
{
   var N = xx_set_detail( geotc, this);
   var S = xx_choose_ring_pattern( N);
   var dz = this.height2 - this.height1;
   var dr = this.radius2 - this.radius1;
   var p = 0;

   for( var n = 0; n < N; n++)
   {
      var x = p + 0.5*S[n];        // Center of the element
      var r = this.radius1 + x * dr;
      var z = this.height1 + x * dz;

      if( this.profile == geotc.PROFILE_LOADED)
         var volts = 0.95 * x + 0.35 * x * x - 0.29 * x * x * x;
      else
      if( this.profile == geotc.PROFILE_BARE)
         var volts = x + 0.65 * x * x - 0.64 * x * x * x;
      else
      if( this.profile == geotc.PROFILE_LINEAR)
         var volts = x
      else
         var volts = 1;           // Defaulting here to PROFILE_DC

      new xx_create_ring( geotc, r, z, S[n]*dr, S[n]*dz, volts);
      p += S[n];
   }

   this.S = S;
}
 
//  Install the secondary.  Will take cylinder, cone, or flat coils.
//
//  The function assumes that the (radius1,height1) end of the coil is
//  to be the grounded end.

function XX_install_secondary( radius1, height1, radius2, height2, 
                               turns, wirad, profile)
{
   var c = this.install_coil( radius1, height1, radius2, height2,
                              turns, wirad);
   c.type = this.TYPE_SEC;
   c.nrings = this.NG_SECONDARY;
   c.profile = profile;
   c.setup = xx_setup_secondary;

   return c;
}

//  Decompose a primary into charge rings.

function xx_setup_primary( geotc)
{
   var N = xx_set_detail( geotc, this);
   var S = xx_choose_ring_pattern( N);
   this.S = S;

   var dr = this.radius2 - this.radius1;
   var dz = this.height2 - this.height1;
   var p = 0;

   for( var n = 0; n < N; n++)
   {
      var x = p + 0.5*S[n];
      var r = this.radius1 + x * dr;
      var z = this.height1 + x * dz;

      new xx_create_ring( geotc, r, z, S[n]*dr, S[n]*dz, 0);
      p += S[n];
   }
}

//  Install a primary.  Cone, cylinder, or flat.  The primary is just
//  another grounded electrode as far as the cap determination is 
//  concerned.

function XX_install_primary( radius1, height1, radius2, height2, turns, tap, wirad)
{
   var c = this.install_coil( radius1, height1, radius2, height2, turns, wirad);

   c.type = this.TYPE_PRI;
   c.nrings = this.NG_PRIMARY;
   c.setup = xx_setup_primary;
   c.tap = tap;

   return c;
}

//  Decompose a toroid into charge rings.

function xx_setup_toroid( geotc)
{
   var N = xx_set_detail( geotc, this);
   var w = (this.outer_radius + this.inner_radius)/2;
   var rd = (this.outer_radius - this.inner_radius)/2;
   var hp = 2 * Math.PI/N;

   var volts = (this.connection == geotc.CONNECT_TOPLOAD) ? 1 : 0;

   for( var i=0; i<N; i++)
   {
      var t = i * hp;
      var hr = rd * hp * Math.sin(t);
      var hz = rd * hp * Math.cos(t);
      if( hr < 0.0) hr = -hr;
      if( hz < 0.0) hz = -hz;
      r = w + rd * Math.cos(t);
      z = this.height + rd * Math.sin(t);
      new xx_create_ring( geotc, r, z, hr, hz, volts);
   }
}

//  Install a toroid.  This can be connected to ground or the secondary
//  topvolts.

function XX_install_toroid( inner_radius, outer_radius, height, connection)
{
   var c = xx_create_component( this, this.TYPE_TOROID,
                                      this.NG_TOROID);
   c.inner_radius = inner_radius;
   c.outer_radius = outer_radius;
   c.height = height;
   c.connection = connection;
   c.setup = xx_setup_toroid;

   return c;
}

//  Decompose a sphere into charge rings.

function xx_setup_sphere( geotc)
{
   var N = xx_set_detail( geotc, this);
   var hp = Math.PI/N;
   var volts = (this.connection == geotc.CONNECT_TOPLOAD) ? 1 : 0;

   for( var i=0; i<N; i++)
   {
      var t = (i+0.5) * hp;
      var hr = hp * this.hradius * Math.cos(t);
      var hz = hp * this.vradius * Math.sin(t);
      if( hr < 0.0) hr = -hr;
      if( hz < 0.0) hz = -hz;
      var r = this.hradius * Math.sin(t);
      var z = this.height + this.vradius * Math.cos(t);
      new xx_create_ring( geotc, r, z, hr, hz, volts);
   }
}

//  Install a sphere topload.  The function is coded for separate
//  horizontal and vertical radii in case someone has a spheroid.
//  connection specifies a grounded sphere or a topvolts sphere.

function XX_install_sphere( hradius, vradius, height, connection)
{
   var c = xx_create_component( this, this.TYPE_SPHERE,
                                      this.NG_SPHERE);
   c.hradius = hradius;
   c.vradius = vradius;
   c.height = height;
   c.connection = connection;
   c.setup = xx_setup_sphere;

   return c;
}

//  Decompose a disc into charge rings.

function xx_setup_disc( geotc)
{
   var N = xx_set_detail( geotc, this); 
   var hr = (this.outer_radius-this.inner_radius)/N;
   var volts = (this.connection == geotc.CONNECT_TOPLOAD) ? 1 : 0;

   if( hr != 0.0)
      for( var k=0; k<N; k++)
         new xx_create_ring( geotc, this.inner_radius + (k+0.5)*hr,
                          this.height, hr, 0.0, volts);
   else
      new xx_create_ring( geotc, this.inner_radius, this.height, 0.0,
                       this.inner_radius * 0.01, volts);
}

//  Install an arbitrary disc shaped electrode.  connection is either topload
//  or ground. To make a strike ring, just call with inner_radius=outer_radius,
//  in which case the program uses a single ring 1cm wide.

function XX_install_disc( inner_radius, outer_radius, height, connection)
{
   var c = xx_create_component( this, this.TYPE_DISC,
                                      this.NG_DISC);
   c.inner_radius = inner_radius;
   c.outer_radius = outer_radius;
   c.height = height;
   c.connection = connection;
   c.setup = xx_setup_disc;

   return c;
}

//  Decompose a cylinder into charge rings.

function xx_setup_cylinder( geotc)
{
   var N = xx_set_detail( geotc, this); 
   var hz = (this.height2-this.height1)/N;
   var volts = (this.connection == geotc.CONNECT_TOPLOAD) ? 1 : 0;

   for( var k=0; k<N; k++) 
      new xx_create_ring( geotc, this.radius, this.height1 + (k+0.5)*hz, 0.0, hz, volts);
}

//  Install an arbitrary cylindrical electrode, eg a core, shield,
//  or a vertical rod, etc.  It is an open cylinder, if you want to
//  plug the ends, use a couple of XX_install_discs().
//
//  As with discs, connection specifies topload or ground.

function XX_install_cylinder( radius, height1, height2, connection)
{
   var c = xx_create_component( this, this.TYPE_CYL,
                                      this.NG_CYL);
   c.radius = radius;
   c.height1 = height1;
   c.height2 = height2;
   c.connection = connection;
   c.setup = xx_setup_cylinder;
   c.ringed = false;
   return c;
}

//  Add a cylindrical, grounded wall around the system.  This is just
//  a wrapper around a grounded cylinder.

function XX_install_wall( radius, height)
{
   var t = this.install_cylinder( radius, 0, height, this.CONNECT_GROUND);
   t.type = this.TYPE_WALL;
   return t;
}

//  Put a grounded roof over the system.   Just a wrapper for a grounded
//  disc.

function XX_install_roof( radius, height)
{
   var t = this.install_disc( 0, radius, height, this.CONNECT_GROUND);
   t.type = this.TYPE_ROOF;
   return t;
}

///////////////////////////////////////////////////////////////////////////
//  Capacitance Stuff.                                                   //
///////////////////////////////////////////////////////////////////////////

function xx_set_detail( geotc, component)
{
   return Math.floor( 0.5 + geotc.detail_val * component.nrings);
}

//  Create a charge ring and add it to the list.

function xx_create_ring( geotc, r, z, hr, hz, volts)
{
   //  Create and initialise a boundary ring and add it to the list.
   geotc.rings[ geotc.nrings++] = this;

   if( hr < 0) hr *= -1;
   if( hz < 0) hz *= -1;
   if( r <= 0) bailout( "faulty ring radius: " + r);

   this.r = r;
   this.z = z;
   this.hr = hr;
   this.hz = hz;
   this.volts = volts;

   this.rc = Math.sqrt(r*r + hr*hr/4);
   if( hr != 0) this.zc =z + (this.rc-r)/hr * hz;
   else this.zc = z;

   //  Calculate the area of each ring, and the electrostatic self potential.

   var e = (this.hr + this.hz)/2;
   if( e <= 0) bailout( "bad ring dimension");
   this.ring_area = 2 * Math.PI * this.r * e;
   this.self_potential = e * Math.log( 1 + 1.414214)/(Math.PI * EPSILON);
}

//  Carlson's elliptic integral - straight out of 'Numerical Recipes':

var ERRTOL = 0.0025

var TINY = 1.5e-38
var BIG = 3.0e37

var C1 = (1.0/24.0)
var C2 = 0.1       
var C3 = (3.0/44.0)
var C4 = (1.0/14.0)

function FMIN( a, b)
{
   if( a < b) return a;
   return b;
}
function FMAX( a, b)
{
   if( a > b) return a;
   return b;
}
function FABS( a)
{
   if( a < 0) return -a;
   return a;
}

function xx_carlson( x, y, z)
{
   var alamb, ave, delx, dely, delz, e2, e3, sqrtx, sqrty, sqrtz, xt, yt, zt;

   if( FMIN( FMIN( x,y), z) < 0.0 || FMIN( FMIN( x+y, x+z), y+z) < TINY ||
       FMAX( FMAX( x,y), z) > BIG) bailout( "invalid arguments in carlson");

   xt = x; yt = y; zt = z; 

   do {
      sqrtx = Math.sqrt( xt); sqrty = Math.sqrt( yt); sqrtz = Math.sqrt( zt);

      alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz;

      xt = 0.25 * (xt+alamb); yt = 0.25 * (yt+alamb); zt = 0.25 * (zt+alamb);

      ave = (xt+yt+zt)/3;

      delx = (ave-xt)/ave; dely = (ave-yt)/ave; delz = (ave-zt)/ave;
   } 
    while( FMAX( FMAX( FABS( delx), FABS( dely)), FABS( delz)) > ERRTOL);

   e2 = delx*dely - delz*delz;
   e3 = delx*dely * delz;

   return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/Math.sqrt(ave);
}

function xx_legendre_elliptic_K( p)
{
   return xx_carlson( 0, 1 - p*p, 1);
}

//  Compute the potential coefficient between two charge rings.

function xx_potential_coefficient_mutual( d, s)
{
   var z = s.z - d.z;
   var t = s.r + d.r;
   var a = Math.sqrt(t*t + z*z);
   var k = 2 * Math.sqrt( s.r * d.r)/a;

   return xx_legendre_elliptic_K( k)/(2*Math.PI*Math.PI*EPSILON*a);
}

//  Compute the potential coefficient between a charge ring and itself.

function xx_potential_coefficient_self( b)
{
   var NSP = 10
   var i;
   var w = ( b.hr*b.hr+b.hz*b.hz)/(NSP*NSP);
   var dr = b.hr/NSP;
   var dz = b.hz/NSP;
   var pc = 0;

   for( i=0; i<NSP; i++)
   {
      var zi = b.z - b.hz/2 + (i+0.5)*dz
      var ri = b.r - b.hr/2 + (i+0.5)*dr;
      var za = zi - b.zc;

      if( za*za + (ri-b.rc)*(ri-b.rc) < w)
      {
         var e = b.rc*b.rc + b.rc * dr + dr*dr/2;
         var z = Math.sqrt(dz*dz + dr*dr)/Math.PI/2;
         var a = Math.sqrt(4*e + z*z);
         var k = 2 * Math.sqrt( e)/a;
         pc += xx_legendre_elliptic_K( k)/a;
      }
      else
      {
         var t = ri + b.rc;
         var a = Math.sqrt(t*t + za*za);
         var k = 2 * Math.sqrt( ri * b.rc)/a;
         pc += xx_legendre_elliptic_K( k)/a;
      }
   }

   return pc/(2*Math.PI*Math.PI*EPSILON*NSP);
}

//  Take the array of charge rings and create a potential matrix.

function xx_setup_potential_matrix( geotc)
{
   if( xx_DEBUG) document.writeln( "computing potential matrix...");

   //  Do the integrations necessary to obtain the coefficients which relate
   //  the potentials at the collocation points to their source charges on the
   //  rings.

   var rings = geotc.rings;
   var N = geotc.nrings;
   var PMAT = geotc.PMAT;

   for( var i=0; i<N; i++)   // Loop over each collocation point...
   {
      if( geotc.progress != null) geotc.progress( Math.floor( 100*i/N));
      if( i >= geotc.pmat_break) PMAT[i+1] = new Array( N+1);
 
      for( var j=0; j<N; j++)       // For each ring of source charge...
      {
         if( i < geotc.pmat_break && j < geotc.pmat_break) continue;

         var H = (rings[j].z - rings[i].z)/rings[i].r;
         var R = rings[j].r/rings[i].r;
         var c;

         if( i == j)
              c = xx_potential_coefficient_self( rings[i]);
         else
              c = xx_potential_coefficient_mutual( rings[i], rings[j]);

         if( (PMAT[i+1][j+1] = c) <= 0.0) 
         {
            bailout( "faulty potential coefficient");
            return false;
         }
      }
   }

   if( xx_DEBUG) document.writeln( "potential matrix setup ok");
   geotc.pmat_break = N;

   return true;
}

//  Multiply a vector by a square matrix.

function xx_real_mulvec(  N, a, d, s)
{
   for( var r=1; r<=N; r++)
   {
      d[r] = 0;
      for( var c=1; c<=N; c++) d[r] += a[r][c] * s[c];
   }
}

//  Return the inner product of two vectors.

function xx_real_inner( N, v1, v2)
{
   var t = 0;

   for( var i=1; i<=N; i++) t += v1[i] * v2[i];
   return t;
}

//  Generalised minimum residual:  Find the vector q which generates the
//  vector v when passed through the matrix pp.  Watch out - the arrays
//  in this function are 1-based.
 
function xx_gmres( pp, q, v, N)
{
   var iter;
   var rnorm;
   var norm;
   var hi; 
   var hip1; 
   var length;

   var ap = new Array( N+1);
   var p = new Array( N+1);

   var c = new Array( N+1);
   var s = new Array( N+1);
   var g = new Array( N+1);
   var y = new Array( N+1);
   var bv = new Array( GMRES_MAXITER+1);
   var bh = new Array( GMRES_MAXITER+1);
 
   rnorm = Math.sqrt( xx_real_inner( N, v, v));

   for( var i=1; i <= N; i++) 
   {
      p[i] = v[i] / rnorm;
      g[i] = 0.0;
   }
   g[1] = rnorm;

   for( var iter = 1; iter <= GMRES_MAXITER && rnorm > GMRES_TOL; iter++)
   {
      bv[iter] = new Array( N+1);
      bh[iter] = new Array( N+1);
      for( var i=1; i <= N; i++) bv[iter][i] = p[i];
    
      xx_real_mulvec( N, pp, ap, p); 

      for( var i=1; i <= N; i++) p[i] = ap[i];
    
      for( var j=1; j <= iter; j++) 
      {
         hi = xx_real_inner( N, ap, bv[j]);
         for( var i=1; i <= N; i++) p[i] -= hi * bv[j][i];
         bh[iter][j] = hi;
      }
    
      norm = Math.sqrt( xx_real_inner( N, p, p));

      for( var i=1; i <= N; i++) p[i] /= norm;
      bh[iter][iter+1] = norm;
    
      for( var i=1; i < iter; i++) 
      {
         hi = bh[iter][i];
         hip1 = bh[iter][i+1];
         bh[iter][i] = c[i] * hi - s[i] * hip1;
         bh[iter][i+1] = c[i] * hip1 + s[i] * hi;
      }
    
      hi = bh[iter][iter];
      hip1 = bh[iter][iter+1];
      length = Math.sqrt(hi * hi + hip1 * hip1);
      c[iter] = hi/length;
      s[iter] = -hip1/length;
    
      bh[iter][iter] = c[iter] * hi - s[iter] * hip1;
      bh[iter][iter+1] = c[iter] * hip1 + s[iter] * hi;
      hi = g[iter];
      g[iter] = c[iter] * hi;
      g[iter+1] = s[iter] * hi;    
    
      rnorm = Math.abs(g[iter+1]);
   }

   iter--;

   for( var i=1; i <= iter; i++) y[i] = g[i];
   for( var i = iter; i > 0; i--) 
   {
      y[i] /= bh[i][i];
      for(j = i-1; j > 0; j--) y[j] -= bh[i][j]*y[i];
   }
   for( var i=1; i <= N; i++) 
   {
      q[i] = 0.0;
      for( var j=1; j <= iter; j++) q[i] += y[j] * bv[j][i];
   }

   if(rnorm > GMRES_TOL) bailout( "gmres: failed to converge");
   if( xx_DEBUG) document.writeln( "gmres iterations: " + iter);
}

//  Delete rings R1 to R2 from the potential matrix and the rings array.
//  Parameters, and rings[] are 0 based, whereas PMAT[][] is 1 based.

function xx_delete_rings( geotc, R1, R2)
{
   var PMAT = geotc.PMAT;
   var r1 = R1+1;
   var r2 = R2+1;
   var n = R2 - R1 + 1;                    // Number of rings to remove

   // Remove rows from PMAT
   for( var j=r1; j<=geotc.nrings-n; j++) PMAT[j] = PMAT[j+n];   

   // Remove columns from PMAT
   for( var j=1; j<=geotc.nrings-n; j++)
      for( var k=r1; k<=geotc.nrings-n; k++) PMAT[j][k] = PMAT[j][k+n];   
   
   // Remove elements from the rings array
   for( var j=R1; j<geotc.nrings-n; j++) geotc.rings[j] = geotc.rings[j+n];   
   geotc.nrings -= n;
   geotc.pmat_break -= n;

   return n;
}

function xx_add_rings( geotc, component)
{
   if( component.ringed) return;

   component.begin = geotc.nrings;
   component.setup( geotc);
   component.end = geotc.nrings;
   component.ringed = true;
}

function xx_remove_rings( geotc, component)
{
   if( !component.ringed) return;

   var n = xx_delete_rings( geotc, component.begin, component.end-1);
   component.ringed = false;

   for( var i=0; i<geotc.count; i++)
   {
      var c = geotc.system[i];
      if( c.ringed && c.begin > component.begin)
      {
         c.begin -= n;
         c.end -= n;
      }
   }

   if( component.type == geotc.TYPE_SEC)
      for( var i=0; i<geotc.count; i++)
      {
         var c = geotc.system[i];
         if( c.type == geotc.TYPE_GND) xx_remove_rings( geotc, c);
      }
}

//
//  Create an array of charge rings from the list of component objects.
//  We deal here with a couple of nasty bits, eg the groundplane needs to
//  know about the secondary radius, if any.
//

function xx_create_charge_rings( geotc)
{
   var sec_radius = null;
   for( var i=0; i<geotc.count; i++)
   {
      var component = geotc.system[i];

      if( component.type == geotc.TYPE_SEC)
      {
         sec_radius = (component.radius1 + component.radius2)/2;
         xx_add_rings( geotc, component);
      } 
      else
      if( component.setup != null &&
          component.type != geotc.TYPE_GND) 
             xx_add_rings( geotc, component);
   }

   for( var i=0; i<geotc.count; i++)
   {
      var component = geotc.system[i];
      if( component.type == geotc.TYPE_GND) 
      {
         component.sec_radius = sec_radius;
         xx_add_rings( geotc, component);
      }
   }

   if( xx_DEBUG) document.writeln( "rings: " + geotc.nrings);
}

//  Return the equivalent shunt capacitance of system, this being the
//  distributed capacitance weighted by the voltage distribution assumed
//  for the system and referred to the topvolts.

function XX_compute_Ces()
{
   var args = XX_compute_Ces.arguments;

   // Convert the component objects to charge rings.
   xx_create_charge_rings( this);

   // Form the potential matrix. Most of the work is in here.
   xx_setup_potential_matrix( this);  

   // Setup a voltage vector.
   var V = new Array( this.nrings+1);            // Voltage array
   if( args.length <= 1)
      for( var i = 0; i<this.nrings; i++) V[i+1] = this.rings[i].volts;
   else
   {
      var source = args[0];
      for( var i = 0; i<this.nrings; i++) 
         if( i >= source.begin && i < source.end) 
         {
            V[i+1] = this.rings[i].volts;
            if( V[i+1] == 0) V[i+1] = 1;
         }
         else V[i+1] = 0;
   }

   // Solve for the potential distribution.
   var Q = new Array( this.nrings+1);            // Charge array
   this.Q = Q;         // Cache the computed charge distribution
   for( var i = 0; i<this.nrings; i++) Q[i+1] = 0;

   // Find the charge distribution Q which generates the given potential
   // distribution V when passed through the potential matrix PMAT.
   xx_gmres( this.PMAT, Q, V, this.nrings);

   //  Sum the charges over the source electrodes.  Since we've
   //  used a nominal 1 volt on the topload, the charge is equal
   //  to the capacitance.  Just which charges we sum as 'source'
   //  and 'destination' of the E flux depends on what arguments
   //  this function has been called with.

   var total_dst_charge = 0.0;    // Total charge on ground rings
   var total_src_charge = 0.0;    // Total charge on source rings

   if( args.length == 0)
   {
      //  Compute the total capacitance of the coil and all the
      //  topload-connected electrodes, to all the grounded electrodes.
      //  Anything with a voltage on it is source, the rest are destination.
   
      for( var k=0; k<this.nrings; k++)
      {
         if( this.rings[k].volts == 0) total_dst_charge += Q[k+1];
         else total_src_charge += Q[k+1];
      }

      if( xx_DEBUG)   
         document.writeln( "source charge: " + total_src_charge*1e12 + 
                           " dest charge: " + total_dst_charge*1e12);

      return total_src_charge;
   }
   else
   if( args.length == 1)
   {
      //  Compute the total capacitance of the single object to 
      //  everything else.  
      var source = args[0];
   
      for( var k=0; k<source.begin; k++) total_dst_charge += Q[k+1];
      for( var k=source.begin; k<source.end; k++) total_src_charge += Q[k+1];
      for( var k=source.end; k<this.nrings; k++) total_dst_charge += Q[k+1];
   
      if( xx_DEBUG)   
         document.writeln( "source charge: " + total_src_charge*1e12 + 
                           " dest charge: " + total_dst_charge*1e12);

      return total_src_charge;
   }
   else
   if( args.length == 2)
   {
      //  Compute the mutual capacitance between the two objects.  The
      //  first component is the source object, at 1 volts, and the
      //  second object is the destination.  This time we are interested
      //  in the destination charge so that we only get the flux arriving
      //  at the second object. If we summed over the source component, we
      //  would get the total C of the source as in the above case.
      var source = args[0];
      var dest = args[1];
   
      for( var k=source.begin; k<source.end; k++) total_src_charge += Q[k+1];
      for( var k=dest.begin; k<dest.end; k++) total_dst_charge += Q[k+1];
   
      if( xx_DEBUG)   
         document.writeln( "source charge: " + total_src_charge*1e12 + 
                           " dest charge: " + total_dst_charge*1e12);

      return -total_dst_charge;
   }
}

///////////////////////////////////////////////////////////////////////////
//  Odds and sods.                                                       //
///////////////////////////////////////////////////////////////////////////

//  Remove a component object from the system.

function XX_remove( c)
{
   xx_remove_rings( this, c);

   for( var i=0; i<this.count; i++)
      if( this.system[i].serial == c.serial)
      {
         for( j=i; j<this.count-1; j++) this.system[j] = this.system[j+1];
         if( --this.count < 0) 
         {
            bailout( "system error (1)");
            return false;
         }
         return true;
      }
 
   if( xx_DEBUG) document.writeln( "component not in system"); 
   return false;
}

//  Empty out the system array.

function XX_reset()
{
   this.count = 0;
   this.timer = (new Date()).getTime();
   this.nrings = 0;
   this.pmat_break = 0;
}

//  Return the elapsed time since the last XX_reset().

function XX_elapsed()
{
   return (new Date()).getTime() - this.timer;
}

// Switch to a new detail setting.  We have to destroy our cache of
// potential coefficients.
 
function XX_detail( d)
{
   oldval = this.detail_val;
   this.detail_val = d; 

   for( var i=0; i<this.count; i++) xx_remove_rings( this, this.system[i]);
   return oldval;
}

//  A funnel for errors. This needs to throw something somewhere but Javascript
//  doesn't seem to have the mechanism.

function bailout( msg)
{
   if( xx_DEBUG) document.writeln( "bailout: " + msg);

   if( !xx_errwin != null)
   {
      if( xx_open_errwin != null) xx_errwin = xx_open_errwin();
   }

   if( xx_errwin != null)
      xx_errwin.document.writeln( "bailout: " + msg + "<br>"); 
   else
   if( !xx_DEBUG)
      document.writeln( "bailout: " + msg);
}

///////////////////////////////////////////////////////////////////////////
//  Compute resonance and VI profiles                                    //
///////////////////////////////////////////////////////////////////////////
//
//

//  Create a 1-based square matrix to represent the coefficients in pn2511
//  eqs 6.3 and 6.4.
// Elements 1..N to represent the voltages at the top of each ring.
// Elements N+1..N+N to represent the currents leaving the top of each ring. 
// Element 2*N+1 to represent the base current.

function xx_enter_Vcoeff( a, n, c)
{
   a[n] += c/2;
   if( n > 1) a[n-1] += c/2;
}

function xx_enter_Icoeff( N, a, n, c)
{
   a[N+n] += c/2;
   if( n > 1) a[N+n-1] += c/2;
   else a[2*N+1] += c/2;
}

function xx_construct_A( N, sec, omega, M, Cint, Cext, Ctor)
{
   var NE = 2*N+1;
   var A = new Array( NE+1);
   for( var i=1; i<=NE; i++) 
   {
      A[i] = new Array( NE+1);
      for( var j=1; j<=NE; j++) A[i][j] = 0;
   }

   // Insert entries for equ 6.3
   for( var n=1; n<=N; n++)
   {
      A[n][n] -= 1;
      if( n > 1) A[n][n-1] += 1;

      for( var m=1; m<=N; m++)
         xx_enter_Icoeff( N, A[n], m, omega * M[n-1][m-1]);
   }

   // Insert entries for equ 6.4.
   for( var n=1; n<=N; n++)
   {
      A[N+n][N+n] -= 1;
      if( n > 1) A[N+n][N+n-1] += 1;
      else A[N+n][2*N+1] += 1;

      xx_enter_Vcoeff( A[N+n], n, -omega * Cext[n-1]);

      for( var m=1; m<=N; m++)
      {
         xx_enter_Vcoeff( A[N+n], n, -omega * Cint[n-1][m-1]);
         xx_enter_Vcoeff( A[N+n], m, omega * Cint[n-1][m-1]);
      }

      xx_enter_Vcoeff( A[N+n], n, -omega * Ctor[n-1]);
      xx_enter_Vcoeff( A[N+n], N, omega * Ctor[n-1]);
   }

   A[2*N+1][2*N+1] = 1;    //  Insert a base current of 1 amp.
  
   return A;
}
 
function XX_compute_F()
{
   var sec = null;
   for( var i = 0; i<this.count; i++)
      if( this.system[i].type == this.TYPE_SEC) sec = this.system[i];

   var result = new Object();

   result.Ldc = this.compute_Ldc( sec);
   sec.profile = this.PROFILE_DC;

   xx_remove_rings( this, sec);
   result.Cdc = this.compute_Ces();
   var N = sec.S.length;
   result.S = sec.S;
   result.Ctop = 0;

   for( var i=0; i<this.nrings; i++) 
      if( (i < sec.begin || i >= sec.end) && 
          this.rings[i].volts != 0) result.Ctop += this.Q[i+1];

   // Create a capacitance matrix.
   var Cext = new Array( N);
   var Cint = new Array( N);
   var Ctor = new Array( N);
   for( var n=0; n<N; n++) Cint[n] = new Array( N);

   for( var i=0; i<N; i++) Cext[i] = this.Q[sec.begin+i+1];

   for( var n=0; n<N; n++)
   { 
      // Voltage array.  Set section 'n' of the secondary to 1 volt, 
      // all the others to zero.
      var V = new Array( this.nrings+1);          
      for( var i = 0; i<this.nrings; i++) V[i+1] = 0;
      V[sec.begin+n+1] = 1;

      // Solve for the resulting charge distribution.
      var Q = new Array( this.nrings+1);            // Charge array
      for( var i = 0; i<this.nrings; i++) Q[i+1] = 0;
      xx_gmres( this.PMAT, Q, V, this.nrings);

      for( var i=0; i<N; i++) Cint[n][i] = -Q[sec.begin+i+1];
      // XXX for( var i=0; i<N; i++) Cint[n][i] = 0;
      Ctor[n] = 0;
      for( var i=0; i<this.nrings; i++) 
         if( (i < sec.begin || i >= sec.end) && 
            this.rings[i].volts != 0) Ctor[n] -= Q[i+1];
   }

   // Create an inductance matrix.
   var M = new Array( N);
   var ND = 3 * this.ND * this.detail_val;
   var H = sec.height2 - sec.height1;
   var R = sec.radius2 - sec.radius1;
   var dH = H/ND;
   var dR = R/ND;

   var nbase = 0;
   for( var n=0; n<N; n++)
   {
      M[n] = new Array( N); 

      var mbase = 0;
      for( var m=0; m<N; m++)
      {
         var nN = Math.floor( 0.5 + ND * sec.S[n]);
         var mN = Math.floor( 0.5 + ND * sec.S[m]);

         var mt = 0;
         for( var j=0; j<nN; j++)
            for( var k=0; k<mN; k++)
            {
               var d = nbase * H + (j + 0.5) * dH - 
                       mbase * H - (k + 0.5) * dH;

               var r1 = (j + 0.5)*dR + sec.radius1 + nbase * R;
               var r2 = (k + 0.5)*dR + sec.radius1 + mbase * R;

               mt += xx_filament( r1, r2, d, sec.wirad);
            }

         M[n][m] = mt * sec.turns * sec.turns/(ND*ND);
         mbase += sec.S[m];
      }

      nbase += sec.S[n];
   }

   var NE = 2*N+1;
   var V = new Array( NE+1);
   var R = new Array( NE+1);

   //  Use Ldc and Cdc to get a start frequency which is guaranteed to be
   //  below Fres.   Hmm, this may not always work, eg with flat spirals, if
   //  Les >> Ldc.  Lets see how it goes.

   var Fstart = 1/(2*Math.PI*Math.sqrt( result.Ldc * result.Cdc));
   var Fstep = 1.2;           // Use a frequency step starting at 20%

   if( xx_DEBUG)
      document.writeln( "Fstart="+Fstart+", initial Cdc="+result.Cdc+
                        ", initial Ldc="+result.Ldc);

   var minp = null
   var minf = 0;
   for( var freq=Fstart; true; freq *= Fstep)
   {
      var omega = 2 * Math.PI * freq;

      var A = xx_construct_A( N, sec, omega, M, Cint, Cext, Ctor);
      for( var n=1; n<=NE; n++) R[n] = 0;
      R[2*N+1] = 1;
      xx_gmres( A, V, R, NE);

      //  Sum the current 'x' which is the total current entering the topload.
      var x = V[2*N] - V[N]*omega*result.Ctop;
      for( var i=0; i<N; i++) x -= (V[N] - V[i+1])*omega*Ctor[i];
      if( xx_DEBUG)
         document.writeln( "freq: "+freq+" "+V[N]+" "+V[2*N]+" x="+x);

      if( (minp == null || V[N] > minp) && x >= 0) 
      { minp = V[N]; minf = freq; }

      //  Break out when the frequency is within 0.1% of the zero crossing.
      //  Otherwise, move to a smaller step.
      if( x < 0) 
      {
         if( Fstep < 1.001) break;     // Finished
         freq /= Fstep; 
         Fstep = Math.sqrt( Fstep);
      }
   }

    //  Our computed Fres.
   freq = minf; 
   var omega = 2 * Math.PI * freq;

   //  Now do a final recompute at this Fres to get the VI profiles.
   var A = xx_construct_A( N, sec, omega, M, Cint, Cext, Ctor);
   for( var n=1; n<=NE; n++) R[n] = 0;
   R[2*N+1] = 1;
   xx_gmres( A, V, R, NE);

   result.N = N;
   result.V = new Array( N);
   result.I = new Array( N);
   result.F = minf;
   for( var n=0; n<N; n++) result.V[n] = V[n+1];
   for( var n=0; n<N; n++) result.I[n] = V[N+n+1];

   result.Cdc = 0;
   for( var n=0; n<N; n++) result.Cdc += Cext[n];
   result.Cdc += result.Ctop;

   //  Calculate Ces and Cee from the voltage profiles. The physical
   //  capacitance is weighted by the average voltage on each element
   //  to obtain the effective capacitance.
   result.Ces = 0;
   for( var n=0; n<N; n++) 
   {
      if( n > 0) result.Ces += (V[n+1]+V[n])/2 * Cext[n];
      else result.Ces += V[n+1]/2 * Cext[n];
   }
   result.Ces += V[N] * result.Ctop;
   result.Ces /= V[N];

   result.Cee = 0;
   for( var n=1; n<=N; n++) 
   {
      if( n > 1) var vn = (V[n]+V[n-1])/2;
      else var vn = V[n]/2;
      result.Cee += Cext[n-1] * vn * vn;

      for( var m=1; m<=N; m++) 
      {
         if( m > 1) var vm = (V[m]+V[m-1])/2;
         else var vm = V[m]/2;
         var v = vn - vm;
         result.Cee += Cint[n-1][m-1] * v * v/2;
      }

      var v = V[n] - V[N];
      result.Cee += Ctor[n-1] * v * v/2;
   }
   result.Cee += V[N]*V[N]*result.Ctop;
   result.Cee /= V[N]*V[N];

   //  Calculate Les and Lee from the current profiles.  The physical
   //  inductance is weighted by the average current in each element
   //  to obtain the effective inductance.
   result.Les = 0;
   for( var n=0; n<N; n++) 
      for( var m=0; m<N; m++) 
      {
         if( m > 0) result.Les += (V[N+m+1]+V[N+m])/2 * M[n][m];
         else result.Les += (V[N+m+1]+V[2*N+1])/2 * M[n][m];
      }

   result.Lee = 0;
   for( var n=1; n<=N; n++) 
      for( var m=1; m<=N; m++) 
      {
         if( m > 1) var Im = (V[N+m]+V[N+m-1])/2;
         else var Im = (V[N+m]+V[2*N+1])/2;
         if( n > 1) var In = (V[N+n]+V[N+n-1])/2;
         else var In = (V[N+n]+V[2*N+1])/2;
         result.Lee += M[n-1][m-1] * Im * In;
      }

   return result;
}

function XX_compute_deck()
{
   var sec = null;
   for( var i = 0; i<this.count; i++)
      if( this.system[i].type == this.TYPE_SEC) sec = this.system[i];

   var result = new Object();

   result.Ldc = this.compute_Ldc( sec);
   sec.profile = this.PROFILE_DC;

   xx_remove_rings( this, sec);
   result.Cdc = this.compute_Ces();
   var N = sec.S.length;
   result.S = sec.S;
   result.Ctop = 0;

   for( var i=0; i<this.nrings; i++) 
      if( (i < sec.begin || i >= sec.end) && 
          this.rings[i].volts != 0) result.Ctop += this.Q[i+1];

   // Create a capacitance matrix.
   var Cext = new Array( N);
   var Cint = new Array( N);
   var Ctor = new Array( N);
   for( var n=0; n<N; n++) Cint[n] = new Array( N);

   for( var i=0; i<N; i++) Cext[i] = this.Q[sec.begin+i+1];

   for( var n=0; n<N; n++)
   { 
      var V = new Array( this.nrings+1);            // Voltage array
      for( var i = 0; i<this.nrings; i++) V[i+1] = 0;
      
      V[sec.begin+n+1] = 1;

      // Solve for the potential distribution.
      var Q = new Array( this.nrings+1);            // Charge array
      for( var i = 0; i<this.nrings; i++) Q[i+1] = 0;
      xx_gmres( this.PMAT, Q, V, this.nrings);

      for( var i=0; i<N; i++) Cint[n][i] = -Q[sec.begin+i+1];
      Ctor[n] = 0;
      for( var i=0; i<this.nrings; i++) 
         if( (i < sec.begin || i >= sec.end) && 
            this.rings[i].volts != 0) Ctor[n] -= Q[i+1];
   }

   // Create an inductance matrix.
   var M = new Array( N);
   var ND = 3 * this.ND * this.detail_val;
   var H = sec.height2 - sec.height1;
   var R = sec.radius2 - sec.radius1;
   var dH = H/ND;
   var dR = R/ND;

   var nbase = 0;
   for( var n=0; n<N; n++)
   {
      M[n] = new Array( N); 

      var mbase = 0;
      for( var m=0; m<N; m++)
      {
         var nN = Math.floor( 0.5 + ND * sec.S[n]);
         var mN = Math.floor( 0.5 + ND * sec.S[m]);

         var mt = 0;
         for( var j=0; j<nN; j++)
            for( var k=0; k<mN; k++)
            {
               var d = nbase * H + (j + 0.5) * dH - 
                       mbase * H - (k + 0.5) * dH;

               var r1 = (j + 0.5)*dR + sec.radius1 + nbase * R;
               var r2 = (k + 0.5)*dR + sec.radius1 + mbase * R;

               mt += xx_filament( r1, r2, d, sec.wirad);
            }

         M[n][m] = mt * sec.turns * sec.turns/(ND*ND);
         mbase += sec.S[m];
      }

      nbase += sec.S[n];
   }

   var NE = 2*N+1;

   document.writeln( ".OPTIONS NOMOD NOPAGE");
   document.writeln( ".AC LIN 500 1K 1000K");
   document.writeln( ".PRINT AC I(VB) V(" + (N+1) + ")");
   document.writeln( "VB 1 0 DC 0 AC 1");
   document.writeln( "");

   document.writeln( ".SUBCKT TC1 N1 N2");
   document.writeln( "RB N1 1 0.1");

   // Inductance matrix...
   for( var n=0; n<N; n++)
   {
      document.writeln( "L" + (n+1) + " " + (n+1) + " " + (n+1) + "R " + M[n][n]);
      document.writeln( "R" + (n+1) + " " + (n+1) + "R " + (n+2) + " 18.3");
   }
   for( var n1=0; n1<N; n1++)
      for( var n2=n1; n2<N; n2++)
   {
      if( n2 != n1)
      {
         var me = M[n1][n2]/Math.sqrt(M[n1][n1] * M[n2][n2]);
         document.writeln( "K" + (n1+1) + "x" + (n2+1) + " L" + (n1+1) + " L" + (n2+1) + " " + me);
      }
   }

   // Cap matrix...
   for( n=0; n<N; n++)
   {
      document.writeln( "CE" + (n+1) + " " + (n+1) + " 0 " + Cext[n]);
      document.writeln( "RCE" + (n+1) + " " + (n+1) + " 0 10000K");
   }
   for( var n1=0; n1<N; n1++)
      for( var n2=n1; n2<N; n2++)
      {
         if( n2 != n1)
            document.writeln( "CI" + (n1+1) + "x" + (n2+1) + " " + (n1+1) + " " + (n2+1) + " " + Cint[n1][n2]);
      }
   for( n=0; n<N; n++)
   {
      document.writeln( "CT" + (n+1) + " " + (n+1) + " " + (N+1) + " " + Ctor[n]);
   }

   document.writeln( "CTOP" + " " + (N+1) + " 0 " + result.Ctop);
   document.writeln( "RT " + (N+1) + " N2 " + "0.1");
   document.writeln( ".ENDS TC1");

   document.writeln( ".END<br>");
   return result;
}

///////////////////////////////////////////////////////////////////////////
//  Interface to Java plotting applet.                                   //
///////////////////////////////////////////////////////////////////////////
//
function XX_graph_V( result)
{
   var d = document.vprof;

   if( d != null && d.init)
   {
      d.init();
      d.wbox_flag = true;
      d.bcolor = document.bgColor;
      d.set_title( "Volts - kV");
      var x = 0;
      d.add_point( 0, 0);
      for( var n=0; n<result.N; n++)
      {
         x += result.S[n];
         d.add_point( x, result.V[n]/1000);
      }
      d.repaint();
   }
   else
      if( xx_DEBUG) document.writeln( "Java not enabled for graph_V");
}

function XX_graph_I( result)
{
   var d = document.iprof;

   if( d != null && d.init)
   {
      d.init();
      d.ecor_flag = true;
      d.wbox_flag = true;
      d.bcolor = document.bgColor;
      d.set_title( "Current - Amps");
      var x = 0;
      d.add_point( 0, 1.0);
      for( var n=0; n<result.N; n++)
      {
         x += result.S[n];
         d.add_point( x, result.I[n]);
      }
      d.repaint();
   }
   else
      if( xx_DEBUG) document.writeln( "Java not enabled for graph_I");
}

function XX_draw()
{
   var d = document.geodraw;

   if( d == null || !d.init)
   {
      if( xx_DEBUG) document.writeln( "Java not enabled for drawing");
      return;
   }

   d.init();
   // d.bcolor = document.bgColor;
   d.bcolor = "#ffffff";
   d.reset();

   for( var i = 0; i<this.count; i++)
   {
      if( this.system[i].type == this.TYPE_TOROID)
      {
         var t = this.system[i];
         d.add_toroid( true, t.inner_radius,
                                            t.outer_radius,
                                            t.height);
      }
      else
      if( this.system[i].type == this.TYPE_SPHERE)
      {
         var t = this.system[i];
         d.add_spheroid( true, t.hradius, t.vradius, t.height);
      }
      else
      if( this.system[i].type == this.TYPE_SEC)
      {
         var t = this.system[i];
         d.add_secondary( true, t.radius1, t.height1,
                                               t.radius2, t.height2);
      }
      else
      if( this.system[i].type == this.TYPE_DISC)
      {
         var t = this.system[i];
         d.add_cone( true, t.inner_radius, t.height,
                                    t.outer_radius, t.height);
      }
      else
      if( this.system[i].type == this.TYPE_ROOF)
      {
         var t = this.system[i];
         d.add_cone( false, t.inner_radius, t.height,
                                    t.outer_radius, t.height);
      }
      else
      if( this.system[i].type == this.TYPE_GND)
      {
         var t = this.system[i];
         d.add_cone( false, 0, 0, t.radius, 0);
      }
      else
      if( this.system[i].type == this.TYPE_CYL)
      {
         var t = this.system[i];
         d.add_cone( true, t.radius, t.height1, t.radius, t.height2);
      }
      else
      if( this.system[i].type == this.TYPE_WALL)
      {
         var t = this.system[i];
         d.add_cone( false, t.radius, t.height1, 
                                           t.radius, t.height2);
      }
      else
      if( this.system[i].type == this.TYPE_PRI)
      {
         var t = this.system[i];
         d.add_cone( true, t.radius1, t.height1,
                                    t.radius2, t.height2);
      }
   }

   d.repaint();
}

///////////////////////////////////////////////////////////////////////////
//  Constructor                                                          //
///////////////////////////////////////////////////////////////////////////
//
//  A constructor to package the whole library.

function GeoTC()
{
   this.version = xx_VERSION;

   // Methods for component installation.
   this.install_secondary = XX_install_secondary;
   this.install_primary = XX_install_primary;
   this.install_toroid = XX_install_toroid;
   this.install_sphere = XX_install_sphere;
   this.install_groundplane = XX_install_groundplane;
   this.install_roof = XX_install_roof;
   this.install_wall = XX_install_wall;
   this.install_disc = XX_install_disc;
   this.install_cylinder = XX_install_cylinder;
   this.install_coil = XX_install_coil;
   this.install_ribbon = XX_install_ribbon;
 
   this.detail = XX_detail;
   this.reset = XX_reset;
   this.elapsed = XX_elapsed;
   this.remove = XX_remove;
   this.progress = null;

   //  Methods for calculations.
   this.compute_Ces = XX_compute_Ces;
   this.compute_Mdc = XX_compute_Mdc;
   this.compute_Ldc = XX_compute_Ldc;
   this.compute_F = XX_compute_F;
   this.compute_deck = XX_compute_deck;

   //  Graphing and drawing.
   this.graph_V = XX_graph_V;
   this.graph_I = XX_graph_I;
   this.draw = XX_draw;

   //  Predefined secondary voltage profile selectors.
   this.PROFILE_LOADED = 1;
   this.PROFILE_BARE = 2;
   this.PROFILE_LINEAR = 3;
   this.PROFILE_DC = 4;

   //  Available connections for components.
   this.CONNECT_TOPLOAD = 1;
   this.CONNECT_GROUND = 2;

   this.system = new Array();                 // Array of system components
   this.count = 0;                            // Count of system components
   this.serial = 1;                      // Serial numbers to ID components 
   this.PMAT = new(Array);                          // The potential matrix
   this.timer = 0;

   this.rings = new Array();                       // Array of charge rings
   this.nrings = 0;                               // Number of charge rings
   this.pmat_break = 0;                       // Max valid offset into PMAT

   // Controls for component rings.
   this.detail_val = 1;
   this.NG_SECONDARY = 15;                 // Number of rings for secondary
   this.NG_GND = 30;                 // Number of rings for the groundplane 
   this.NG_DISC = 15;                // Number of rings for disc structures
   this.NG_PRIMARY = 8;                      // Number of rings for primary
   this.NG_TOROID = 15;                     // Number of rings for a toroid 
   this.NG_SPHERE = 15;                     // Number of rings for a sphere
   this.NG_CYL = 15;             // Number of rings for cylinder structures
   this.ND = 150;           // Number of current elements to split windings

   //  Component types.
   this.TYPE_PRI = 1;
   this.TYPE_SEC = 2;
   this.TYPE_SPHERE = 3;
   this.TYPE_TOROID = 4;
   this.TYPE_COIL = 5;
   this.TYPE_DISC = 6;
   this.TYPE_CYL = 7;
   this.TYPE_GND = 8;
   this.TYPE_WALL = 9;
   this.TYPE_ROOF = 10;

   //  Coil types.
   this.COIL_THIN = 1;
   this.COIL_RIBBON = 2;

   //  Setup for error messages if the user has supplied a function to
   //  open an error window.
   var args = GeoTC.arguments;
   if( args.length == 1) xx_open_errwin = args[0];
}

//
//
///////////////////////////////////////////////////////////////////////////
//  End of library functions.                                            //
///////////////////////////////////////////////////////////////////////////
//