File indexing completed on 2024-05-05 03:41:52
0001 /* cktsim.js 0002 * 0003 * SPDX-FileCopyrightText: 2011 Massachusetts Institute of Technology 0004 * SPDX-License-Identifier: AGPL-3.0-only 0005 * 0006 * Authors: 0007 * Massachusetts Institute of Technology (original engine) 0008 * Aiswarya Kaitheri Kandoth <aiswaryakk29@gmail.com> (minimal edits to integrate in GCompris) 0009 * 0010 * Copied from https://github.com/edx/edx-platform/blob/master/common/lib/xmodule/xmodule/js/src/capa/schematic.js 0011 */ 0012 0013 /* eslint-disable */ 0014 0015 0016 ////////////////////////////////////////////////////////////////////////////// 0017 // 0018 // Circuit simulator 0019 // 0020 ////////////////////////////////////////////////////////////////////////////// 0021 0022 // create a circuit for simulation using "new cktsim.Circuit()" 0023 0024 // for modified nodal analysis (MNA) stamps see 0025 // http://www.analog-electronics.eu/analog-electronics/modified-nodal-analysis/modified-nodal-analysis.xhtml 0026 0027 var cktsim = (function() { 0028 /////////////////////////////////////////////////////////////////////////////// 0029 // 0030 // Circuit 0031 // 0032 ////////////////////////////////////////////////////////////////////////////// 0033 0034 // types of "nodes" in the linear system 0035 var T_VOLTAGE = 0; 0036 var T_CURRENT = 1; 0037 0038 var v_newt_lim = 0.3; // Voltage limited Newton great for Mos/diodes 0039 var v_abstol = 1e-6; // Absolute voltage error tolerance 0040 var i_abstol = 1e-12; // Absolute current error tolerance 0041 var eps = 1.0e-12; // A very small number compared to one. 0042 var dc_max_iters = 1000; // max iterations before giving up 0043 var max_tran_iters = 20; // max iterations before giving up 0044 var time_step_increase_factor = 2.0; // How much can lte let timestep grow. 0045 var lte_step_decrease_factor = 8; // Limit lte one-iter timestep shrink. 0046 var nr_step_decrease_factor = 4; // Newton failure timestep shrink. 0047 var reltol = 0.0001; // Relative tol to max observed value 0048 var lterel = 10; // LTE/Newton tolerance ratio (> 10!) 0049 var res_check_abs = Math.sqrt(i_abstol); // Loose Newton residue check 0050 var res_check_rel = Math.sqrt(reltol); // Loose Newton residue check 0051 0052 function Circuit() { 0053 this.node_map = []; 0054 this.ntypes = []; 0055 this.initial_conditions = []; 0056 this.devices = []; 0057 this.device_map = []; 0058 this.voltage_sources = []; 0059 this.current_sources = []; 0060 this.finalized = false; 0061 this.diddc = false; 0062 this.node_index = -1; 0063 this.periods = 1; 0064 // Added for GCompris to store current results 0065 this.GCCurrentResults = []; 0066 this.GCWarning = ""; 0067 } 0068 0069 // alert managed in GCompris 0070 // (note: all alert function calls have been changed to this.alert... + added qsTr() for translation) 0071 Circuit.prototype.alert = function(message_) { 0072 if(this.GCWarning === "") { 0073 this.GCWarning = message_; 0074 } 0075 else { 0076 this.GCWarning = this.GCWarning + "\n\n" + message_; 0077 } 0078 return; 0079 } 0080 0081 // index of ground node 0082 Circuit.prototype.gnd_node = function() { 0083 return -1; 0084 } 0085 0086 // allocate a new node index 0087 Circuit.prototype.node = function(name,ntype,ic) { 0088 this.node_index += 1; 0089 if (name) this.node_map[name] = this.node_index; 0090 this.ntypes.push(ntype); 0091 this.initial_conditions.push(ic); 0092 return this.node_index; 0093 } 0094 0095 // call to finalize the circuit in preparation for simulation 0096 Circuit.prototype.finalize = function() { 0097 if (!this.finalized) { 0098 this.finalized = true; 0099 this.N = this.node_index + 1; // number of nodes 0100 0101 // give each device a chance to finalize itself 0102 for (var i = this.devices.length - 1; i >= 0; --i) 0103 this.devices[i].finalize(this); 0104 0105 // set up augmented matrix and various temp vectors 0106 this.matrix = mat_make(this.N, this.N+1); 0107 this.Gl = mat_make(this.N, this.N); // Matrix for linear conductances 0108 this.G = mat_make(this.N, this.N); // Complete conductance matrix 0109 this.C = mat_make(this.N, this.N); // Matrix for linear L's and C's 0110 0111 this.soln_max = new Array(this.N); // max abs value seen for each unknown 0112 this.abstol = new Array(this.N); 0113 this.solution = new Array(this.N); 0114 this.rhs = new Array(this.N); 0115 for (var i = this.N - 1; i >= 0; --i) { 0116 this.soln_max[i] = 0.0; 0117 this.abstol[i] = this.ntypes[i] == T_VOLTAGE ? v_abstol : i_abstol; 0118 this.solution[i] = 0.0; 0119 this.rhs[i] = 0.0; 0120 } 0121 0122 // Load up the linear elements once and for all 0123 for (var i = this.devices.length - 1; i >= 0; --i) { 0124 this.devices[i].load_linear(this) 0125 } 0126 0127 // Check for voltage source loops. 0128 var n_vsrc = this.voltage_sources.length; 0129 if (n_vsrc > 0) { // At least one voltage source 0130 var GV = mat_make(n_vsrc, this.N); // Loop check 0131 for (var i = n_vsrc - 1; i >= 0; --i) { 0132 var branch = this.voltage_sources[i].branch; 0133 for (var j = this.N - 1; j >= 0; j--) 0134 GV[i][j] = this.Gl[branch][j]; 0135 } 0136 var rGV = mat_rank(GV); 0137 if (rGV < n_vsrc) { 0138 this.alert(qsTr('Warning! Circuit has a voltage source loop or a source shorted by a wire, please remove the wire causing the short.')); 0139 this.alert(qsTr('Warning! Simulator might produce meaningless results or no result with this circuit.')); 0140 return false; 0141 } 0142 } 0143 } 0144 return true; 0145 } 0146 0147 // load circuit from JSON netlist (see schematic.js) 0148 Circuit.prototype.load_netlist = function(netlist) { 0149 // set up mapping for all ground connections 0150 for (var i = netlist.length - 1; i >= 0; --i) { 0151 var component = netlist[i]; 0152 var type = component[0]; 0153 if (type == 'g') { 0154 var connections = component[3]; 0155 this.node_map[connections[0]] = this.gnd_node(); 0156 } 0157 } 0158 0159 // process each component in the JSON netlist (see schematic.js for format) 0160 var found_ground = false; 0161 for (var i = netlist.length - 1; i >= 0; --i) { 0162 var component = netlist[i]; 0163 var type = component[0]; 0164 0165 // ignore wires, ground connections, scope probes and view info 0166 if (type == 'view' || type == 'w' || type == 'g' || type == 's' || type == 'L') { 0167 continue; 0168 } 0169 0170 var properties = component[2]; 0171 var name = properties['name']; 0172 if (name==undefined || name=='') 0173 name = '_' + properties['_json_'].toString(); 0174 0175 // convert node names to circuit indicies 0176 var connections = component[3]; 0177 for (var j = connections.length - 1; j >= 0; --j) { 0178 var node = connections[j]; 0179 var index = this.node_map[node]; 0180 if (index == undefined) index = this.node(node,T_VOLTAGE); 0181 else if (index == this.gnd_node()) found_ground = true; 0182 connections[j] = index; 0183 } 0184 0185 // process the component 0186 if (type == 'r') // resistor 0187 this.r(connections[0],connections[1],properties['r'],name); 0188 else if (type == 'd') // diode 0189 this.d(connections[0],connections[1],properties['area'],properties['type'],name); 0190 else if (type == 'c') // capacitor 0191 this.c(connections[0],connections[1],properties['c'],name); 0192 else if (type == 'l') // inductor 0193 this.l(connections[0],connections[1],properties['l'],name); 0194 else if (type == 'v') // voltage source 0195 this.v(connections[0],connections[1],properties['value'],name); 0196 else if (type == 'i') // current source 0197 this.i(connections[0],connections[1],properties['value'],name); 0198 else if (type == 'o') // op amp 0199 this.opamp(connections[0],connections[1],connections[2],connections[3],properties['A'],name); 0200 else if (type == 'n') // n fet 0201 this.n(connections[0],connections[1],connections[2],properties['W/L'],name); 0202 else if (type == 'p') // p fet 0203 this.p(connections[0],connections[1],connections[2],properties['W/L'],name); 0204 else if (type == 'a') // current probe == 0-volt voltage source 0205 this.v(connections[0],connections[1],'0',name); 0206 } 0207 0208 // if (!found_ground) { // No ground on schematic 0209 // alert('Please make at least one connection to ground (inverted T symbol)'); 0210 // return false; 0211 // } 0212 return true; 0213 } 0214 0215 // if converges: updates this.solution, this.soln_max, returns iter count 0216 // otherwise: return undefined and set this.problem_node 0217 // Load should compute -f and df/dx (note the sign pattern!) 0218 Circuit.prototype.find_solution = function(load,maxiters) { 0219 var soln = this.solution; 0220 var rhs = this.rhs; 0221 var d_sol = []; 0222 var abssum_compare; 0223 var converged,abssum_old=0, abssum_rhs; 0224 var use_limiting = false; 0225 var down_count = 0; 0226 var thresh; 0227 0228 // iteratively solve until values convere or iteration limit exceeded 0229 for (var iter = 0; iter < maxiters; iter++) { 0230 // set up equations 0231 load(this,soln,rhs); 0232 0233 // Compute norm of rhs, assume variables of v type go with eqns of i type 0234 abssum_rhs = 0; 0235 for (var i = this.N - 1; i >= 0; --i) 0236 if (this.ntypes[i] == T_VOLTAGE) 0237 abssum_rhs += Math.abs(rhs[i]); 0238 0239 if ((iter > 0) && (use_limiting == false) && (abssum_old < abssum_rhs)) { 0240 // Old rhsnorm was better, undo last iter and turn on limiting 0241 for (var i = this.N - 1; i >= 0; --i) 0242 soln[i] -= d_sol[i]; 0243 iter -= 1; 0244 use_limiting = true; 0245 } 0246 else { // Compute the Newton delta 0247 d_sol = mat_solve_rq(this.matrix,rhs); 0248 0249 // If norm going down for ten iters, stop limiting 0250 if (abssum_rhs < abssum_old) 0251 down_count += 1; 0252 else 0253 down_count = 0; 0254 if (down_count > 10) { 0255 use_limiting = false; 0256 down_count = 0; 0257 } 0258 0259 // Update norm of rhs 0260 abssum_old = abssum_rhs; 0261 } 0262 0263 // Update the worst case abssum for comparison. 0264 if ((iter == 0) || (abssum_rhs > abssum_compare)) 0265 abssum_compare = abssum_rhs; 0266 0267 // Check residue convergence, but loosely, and give up 0268 // on last iteration 0269 if ( (iter < (maxiters - 1)) && 0270 (abssum_rhs > (res_check_abs+res_check_rel*abssum_compare))) 0271 converged = false; 0272 else converged = true; 0273 0274 0275 // Update solution and check delta convergence 0276 for (var i = this.N - 1; i >= 0; --i) { 0277 // Simple voltage step limiting to encourage Newton convergence 0278 if (use_limiting) { 0279 if (this.ntypes[i] == T_VOLTAGE) { 0280 d_sol[i] = (d_sol[i] > v_newt_lim) ? v_newt_lim : d_sol[i]; 0281 d_sol[i] = (d_sol[i] < -v_newt_lim) ? -v_newt_lim : d_sol[i]; 0282 } 0283 } 0284 soln[i] += d_sol[i]; 0285 thresh = this.abstol[i] + reltol*this.soln_max[i]; 0286 if (Math.abs(d_sol[i]) > thresh) { 0287 converged = false; 0288 this.problem_node = i; 0289 } 0290 } 0291 0292 if (converged == true) { 0293 for (var i = this.N - 1; i >= 0; --i) 0294 if (Math.abs(soln[i]) > this.soln_max[i]) 0295 this.soln_max[i] = Math.abs(soln[i]); 0296 return iter+1; 0297 } 0298 } 0299 return undefined; 0300 } 0301 0302 // DC analysis 0303 Circuit.prototype.dc = function() { 0304 // Allocation matrices for linear part, etc. 0305 if (this.finalize() == false) 0306 return undefined; 0307 0308 // Define -f and df/dx for Newton solver 0309 function load_dc(ckt,soln,rhs) { 0310 // rhs is initialized to -Gl * soln 0311 mat_v_mult(ckt.Gl, soln, rhs, -1.0); 0312 // G matrix is initialized with linear Gl 0313 mat_copy(ckt.Gl,ckt.G); 0314 // Now load up the nonlinear parts of rhs and G 0315 for (var i = ckt.devices.length - 1; i >= 0; --i) 0316 ckt.devices[i].load_dc(ckt,soln,rhs); 0317 // G matrix is copied in to the system matrix 0318 mat_copy(ckt.G,ckt.matrix); 0319 } 0320 0321 // find the operating point 0322 var iterations = this.find_solution(load_dc,dc_max_iters); 0323 0324 if (typeof iterations == 'undefined') { 0325 // too many iterations 0326 if (this.current_sources.length > 0) { 0327 this.alert('Newton Method Failed, do your current sources have a conductive path to ground?'); 0328 } else { 0329 this.alert('Newton Method Failed, it may be your circuit or it may be our simulator.'); 0330 } 0331 0332 return undefined 0333 } else { 0334 // Note that a dc solution was computed 0335 this.diddc = true; 0336 // create solution dictionary 0337 var result = []; 0338 // capture node voltages 0339 for (var name in this.node_map) { 0340 var index = this.node_map[name]; 0341 result[name] = (index == -1) ? 0 : this.solution[index]; 0342 } 0343 // reinit results list for GCompris 0344 this.GCCurrentResults = []; 0345 // capture branch currents from voltage sources 0346 for (var i = this.voltage_sources.length - 1; i >= 0; --i) { 0347 var v = this.voltage_sources[i]; 0348 result['I('+v.name+')'] = this.solution[v.branch]; 0349 this.GCCurrentResults.push(result['I('+v.name+')']); 0350 } 0351 return result; 0352 } 0353 } 0354 0355 // Transient analysis (needs work!) 0356 Circuit.prototype.tran = function(ntpts, tstart, tstop, probenames, no_dc) { 0357 0358 // Define -f and df/dx for Newton solver 0359 function load_tran(ckt,soln,rhs) { 0360 // Crnt is initialized to -Gl * soln 0361 mat_v_mult(ckt.Gl, soln, ckt.c,-1.0); 0362 // G matrix is initialized with linear Gl 0363 mat_copy(ckt.Gl,ckt.G); 0364 // Now load up the nonlinear parts of crnt and G 0365 for (var i = ckt.devices.length - 1; i >= 0; --i) 0366 ckt.devices[i].load_tran(ckt,soln,ckt.c,ckt.time); 0367 // Exploit the fact that storage elements are linear 0368 mat_v_mult(ckt.C, soln, ckt.q, 1.0); 0369 // -rhs = c - dqdt 0370 for (var i = ckt.N-1; i >= 0; --i) { 0371 var dqdt = ckt.alpha0*ckt.q[i] + ckt.alpha1*ckt.oldq[i] + 0372 ckt.alpha2*ckt.old2q[i]; 0373 rhs[i] = ckt.beta0[i]*ckt.c[i] + ckt.beta1[i]*ckt.oldc[i] - dqdt; 0374 } 0375 // matrix = beta0*G + alpha0*C. 0376 mat_scale_add(ckt.G,ckt.C,ckt.beta0,ckt.alpha0,ckt.matrix); 0377 } 0378 0379 var p = new Array(3); 0380 function interp_coeffs(t, t0, t1, t2) { 0381 // Poly coefficients 0382 var dtt0 = (t - t0); 0383 var dtt1 = (t - t1); 0384 var dtt2 = (t - t2); 0385 var dt0dt1 = (t0 - t1); 0386 var dt0dt2 = (t0 - t2); 0387 var dt1dt2 = (t1 - t2); 0388 p[0] = (dtt1*dtt2)/(dt0dt1 * dt0dt2); 0389 p[1] = (dtt0*dtt2)/(-dt0dt1 * dt1dt2); 0390 p[2] = (dtt0*dtt1)/(dt0dt2 * dt1dt2); 0391 return p; 0392 } 0393 0394 function pick_step(ckt, step_index) { 0395 var min_shrink_factor = 1.0/lte_step_decrease_factor; 0396 var max_growth_factor = time_step_increase_factor; 0397 var N = ckt.N; 0398 var p = interp_coeffs(ckt.time, ckt.oldt, ckt.old2t, ckt.old3t); 0399 var trapcoeff = 0.5*(ckt.time - ckt.oldt)/(ckt.time - ckt.old3t); 0400 var maxlteratio = 0.0; 0401 for (var i = ckt.N-1; i >= 0; --i) { 0402 if (ckt.ltecheck[i]) { // Check lte on variable 0403 var pred = p[0]*ckt.oldsol[i] + p[1]*ckt.old2sol[i] + p[2]*ckt.old3sol[i]; 0404 var lte = Math.abs((ckt.solution[i] - pred))*trapcoeff; 0405 var lteratio = lte/(lterel*(ckt.abstol[i] + reltol*ckt.soln_max[i])); 0406 maxlteratio = Math.max(maxlteratio, lteratio); 0407 } 0408 } 0409 var new_step; 0410 var lte_step_ratio = 1.0/Math.pow(maxlteratio,1/3); // Cube root because trap 0411 if (lte_step_ratio < 1.0) { // Shrink the timestep to make lte 0412 lte_step_ratio = Math.max(lte_step_ratio,min_shrink_factor); 0413 new_step = (ckt.time - ckt.oldt)*0.75*lte_step_ratio; 0414 new_step = Math.max(new_step, ckt.min_step); 0415 } else { 0416 lte_step_ratio = Math.min(lte_step_ratio, max_growth_factor); 0417 if (lte_step_ratio > 1.2) /* Increase timestep due to lte. */ 0418 new_step = (ckt.time - ckt.oldt) * lte_step_ratio / 1.2; 0419 else 0420 new_step = (ckt.time - ckt.oldt); 0421 new_step = Math.min(new_step, ckt.max_step); 0422 } 0423 return new_step; 0424 } 0425 0426 // Standard to do a dc analysis before transient 0427 // Otherwise, do the setup also done in dc. 0428 no_dc = false; 0429 if ((this.diddc == false) && (no_dc == false)) { 0430 if (this.dc() == undefined) { // DC failed, realloc mats and vects. 0431 this.alert('DC failed, trying transient analysis from zero.'); 0432 this.finalized = false; // Reset the finalization. 0433 if (this.finalize() == false) 0434 return undefined; 0435 } 0436 } 0437 else { 0438 if (this.finalize() == false) // Allocate matrices and vectors. 0439 return undefined; 0440 } 0441 0442 // Tired of typing this, and using "with" generates hate mail. 0443 var N = this.N; 0444 0445 // build array to hold list of results for each variable 0446 // last entry is for timepoints. 0447 var response = new Array(N + 1); 0448 for (var i = N; i >= 0; --i) response[i] = []; 0449 0450 // Allocate back vectors for up to a second order method 0451 this.old3sol = new Array(this.N); 0452 this.old3q = new Array(this.N); 0453 this.old2sol = new Array(this.N); 0454 this.old2q = new Array(this.N); 0455 this.oldsol = new Array(this.N); 0456 this.oldq = new Array(this.N); 0457 this.q = new Array(this.N); 0458 this.oldc = new Array(this.N); 0459 this.c = new Array(this.N); 0460 this.alpha0 = 1.0; 0461 this.alpha1 = 0.0; 0462 this.alpha2 = 0.0; 0463 this.beta0 = new Array(this.N); 0464 this.beta1 = new Array(this.N); 0465 0466 // Mark a set of algebraic variable (don't miss hidden ones!). 0467 this.ar = this.algebraic(this.C); 0468 0469 // Non-algebraic variables and probe variables get lte 0470 this.ltecheck = new Array(this.N); 0471 for (var i = N; i >= 0; --i) 0472 this.ltecheck[i] = (this.ar[i] == 0); 0473 0474 for (var name in this.node_map) { 0475 var index = this.node_map[name]; 0476 for (var i = probenames.length; i >= 0; --i) { 0477 if (name == probenames[i]) { 0478 this.ltecheck[index] = true; 0479 break; 0480 } 0481 } 0482 } 0483 0484 // Check for periodic sources 0485 var period = tstop - tstart; 0486 for (var i = this.voltage_sources.length - 1; i >= 0; --i) { 0487 var per = this.voltage_sources[i].src.period; 0488 if (per > 0) 0489 period = Math.min(period, per); 0490 } 0491 for (var i = this.current_sources.length - 1; i >= 0; --i) { 0492 var per = this.current_sources[i].src.period; 0493 if (per > 0) 0494 period = Math.min(period, per); 0495 } 0496 this.periods = Math.ceil((tstop - tstart)/period); 0497 0498 this.time = tstart; 0499 // ntpts adjusted by numbers of periods in input 0500 this.max_step = (tstop - tstart)/(this.periods*ntpts); 0501 this.min_step = this.max_step/1e8; 0502 var new_step = this.max_step/1e6; 0503 this.oldt = this.time - new_step; 0504 0505 // Initialize old crnts, charges, and solutions. 0506 load_tran(this,this.solution,this.rhs) 0507 for (var i = N-1; i >= 0; --i) { 0508 this.old3sol[i] = this.solution[i]; 0509 this.old2sol[i] = this.solution[i]; 0510 this.oldsol[i] = this.solution[i]; 0511 this.old3q[i] = this.q[i]; 0512 this.old2q[i] = this.q[i]; 0513 this.oldq[i] = this.q[i]; 0514 this.oldc[i] = this.c[i]; 0515 } 0516 0517 var beta0,beta1; 0518 // Start with two pseudo-Euler steps, maximum 50000 steps/period 0519 var max_nsteps = this.periods*50000; 0520 for(var step_index = -3; step_index < max_nsteps; step_index++) { 0521 // Save the just computed solution, and move back q and c. 0522 for (var i = this.N - 1; i >= 0; --i) { 0523 if (step_index >= 0) 0524 response[i].push(this.solution[i]); 0525 this.oldc[i] = this.c[i]; 0526 this.old3sol[i] = this.old2sol[i]; 0527 this.old2sol[i] = this.oldsol[i]; 0528 this.oldsol[i] = this.solution[i]; 0529 this.old3q[i] = this.oldq[i]; 0530 this.old2q[i] = this.oldq[i]; 0531 this.oldq[i] = this.q[i]; 0532 } 0533 0534 if (step_index < 0) { // Take a prestep using BE 0535 this.old3t = this.old2t - (this.oldt-this.old2t) 0536 this.old2t = this.oldt - (tstart-this.oldt) 0537 this.oldt = tstart - (this.time - this.oldt); 0538 this.time = tstart; 0539 beta0 = 1.0; 0540 beta1 = 0.0; 0541 } else { // Take a regular step 0542 // Save the time, and rotate time wheel 0543 response[this.N].push(this.time); 0544 this.old3t = this.old2t; 0545 this.old2t = this.oldt; 0546 this.oldt = this.time; 0547 // Make sure we come smoothly in to the interval end. 0548 if (this.time >= tstop) break; // We're done. 0549 else if(this.time + new_step > tstop) 0550 this.time = tstop; 0551 else if(this.time + 1.5*new_step > tstop) 0552 this.time += (2/3)*(tstop - this.time); 0553 else 0554 this.time += new_step; 0555 0556 // Use trap (average old and new crnts. 0557 beta0 = 0.5; 0558 beta1 = 0.5; 0559 } 0560 0561 // For trap rule, turn off current avging for algebraic eqns 0562 for (var i = this.N - 1; i >= 0; --i) { 0563 this.beta0[i] = beta0 + this.ar[i]*beta1; 0564 this.beta1[i] = (1.0 - this.ar[i])*beta1; 0565 } 0566 0567 // Loop to find NR converging timestep with okay LTE 0568 while (true) { 0569 // Set the timestep coefficients (alpha2 is for bdf2). 0570 this.alpha0 = 1.0/(this.time - this.oldt); 0571 this.alpha1 = -this.alpha0; 0572 this.alpha2 = 0; 0573 0574 // If timestep is 1/10,000th of tstop, just use BE. 0575 if ((this.time-this.oldt) < 1.0e-4*tstop) { 0576 for (var i = this.N - 1; i >= 0; --i) { 0577 this.beta0[i] = 1.0; 0578 this.beta1[i] = 0.0; 0579 } 0580 } 0581 // Use Newton to compute the solution. 0582 var iterations = this.find_solution(load_tran,max_tran_iters); 0583 0584 // If NR succeeds and stepsize is at min, accept and newstep=maxgrowth*minstep. 0585 // Else if Newton Fails, shrink step by a factor and try again 0586 // Else LTE picks new step, if bigger accept current step and go on. 0587 if ((iterations != undefined) && 0588 (step_index <= 0 || (this.time-this.oldt) < (1+reltol)*this.min_step)) { 0589 if (step_index > 0) new_step = time_step_increase_factor*this.min_step; 0590 break; 0591 } else if (iterations == undefined) { // NR nonconvergence, shrink by factor 0592 this.time = this.oldt + 0593 (this.time - this.oldt)/nr_step_decrease_factor; 0594 } else { // Check the LTE and shrink step if needed. 0595 new_step = pick_step(this, step_index); 0596 if (new_step < (1.0 - reltol)*(this.time - this.oldt)) { 0597 this.time = this.oldt + new_step; // Try again 0598 } 0599 else 0600 break; // LTE okay, new_step for next step 0601 } 0602 } 0603 } 0604 0605 // create solution dictionary 0606 var result = []; 0607 for (var name in this.node_map) { 0608 var index = this.node_map[name]; 0609 result[name] = (index == -1) ? 0 : response[index]; 0610 } 0611 // capture branch currents from voltage sources 0612 for (var i = this.voltage_sources.length - 1; i >= 0; --i) { 0613 var v = this.voltage_sources[i]; 0614 result['I('+v.name+')'] = response[v.branch]; 0615 } 0616 0617 result['_time_'] = response[this.N]; 0618 return result; 0619 } 0620 0621 // AC analysis: npts/decade for freqs in range [fstart,fstop] 0622 // result['_frequencies_'] = vector of log10(sample freqs) 0623 // result['xxx'] = vector of dB(response for node xxx) 0624 // NOTE: Normalization removed in schematic.js, jkw. 0625 Circuit.prototype.ac = function(npts,fstart,fstop,source_name) { 0626 0627 if (this.dc() == undefined) { // DC failed, realloc mats and vects. 0628 return undefined; 0629 } 0630 0631 var N = this.N; 0632 var G = this.G; 0633 var C = this.C; 0634 0635 // Complex numbers, we're going to need a bigger boat 0636 var matrixac = mat_make(2*N, (2*N)+1); 0637 0638 // Get the source used for ac 0639 if (this.device_map[source_name] === undefined) { 0640 this.alert('AC analysis refers to unknown source ' + source_name); 0641 return 'AC analysis failed, unknown source'; 0642 } 0643 this.device_map[source_name].load_ac(this,this.rhs); 0644 0645 // build array to hold list of magnitude and phases for each node 0646 // last entry is for frequency values 0647 var response = new Array(2*N + 1); 0648 for (var i = 2*N; i >= 0; --i) response[i] = []; 0649 0650 // multiplicative frequency increase between freq points 0651 var delta_f = Math.exp(Math.LN10/npts); 0652 0653 var phase_offset = new Array(N); 0654 for (var i = N-1; i >= 0; --i) phase_offset[i] = 0; 0655 0656 var f = fstart; 0657 fstop *= 1.0001; // capture that last freq point! 0658 while (f <= fstop) { 0659 var omega = 2 * Math.PI * f; 0660 response[2*N].push(f); // 2*N for magnitude and phase 0661 0662 // Find complex x+jy that sats Gx-omega*Cy=rhs; omega*Cx+Gy=0 0663 // Note: solac[0:N-1]=x, solac[N:2N-1]=y 0664 for (var i = N-1; i >= 0; --i) { 0665 // First the rhs, replicated for real and imaginary 0666 matrixac[i][2*N] = this.rhs[i]; 0667 matrixac[i+N][2*N] = 0; 0668 0669 for (var j = N-1; j >= 0; --j) { 0670 matrixac[i][j] = G[i][j]; 0671 matrixac[i+N][j+N] = G[i][j]; 0672 matrixac[i][j+N] = -omega*C[i][j]; 0673 matrixac[i+N][j] = omega*C[i][j]; 0674 } 0675 } 0676 0677 // Compute the small signal response 0678 var solac = mat_solve(matrixac); 0679 0680 // Save magnitude and phase 0681 for (var i = N - 1; i >= 0; --i) { 0682 var mag = Math.sqrt(solac[i]*solac[i] + solac[i+N]*solac[i+N]); 0683 response[i].push(mag); 0684 0685 // Avoid wrapping phase, add or sub 180 for each jump 0686 var phase = 180*(Math.atan2(solac[i+N],solac[i])/Math.PI); 0687 var phasei = response[i+N]; 0688 var L = phasei.length; 0689 // Look for a one-step jump greater than 90 degrees 0690 if (L > 1) { 0691 var phase_jump = phase + phase_offset[i] - phasei[L-1]; 0692 if (phase_jump > 90) { 0693 phase_offset[i] -= 360; 0694 } else if (phase_jump < -90) { 0695 phase_offset[i] += 360; 0696 } 0697 } 0698 response[i+N].push(phase + phase_offset[i]); 0699 } 0700 f *= delta_f; // increment frequency 0701 } 0702 0703 // create solution dictionary 0704 var result = []; 0705 for (var name in this.node_map) { 0706 var index = this.node_map[name]; 0707 result[name] = (index == -1) ? 0 : response[index]; 0708 result[name+'_phase'] = (index == -1) ? 0 : response[index+N]; 0709 } 0710 result['_frequencies_'] = response[2*N]; 0711 return result; 0712 } 0713 0714 // Helper for adding devices to a circuit, warns on duplicate device names. 0715 Circuit.prototype.add_device = function(d,name) { 0716 // Add device to list of devices and to device map 0717 this.devices.push(d); 0718 d.name = name; 0719 if (name) { 0720 if (this.device_map[name] === undefined) 0721 this.device_map[name] = d; 0722 else { 0723 this.alert('Warning: two circuit elements share the same name ' + name); 0724 this.device_map[name] = d; 0725 } 0726 } 0727 return d; 0728 } 0729 0730 Circuit.prototype.r = function(n1,n2,v,name) { 0731 // try to convert string value into numeric value, barf if we can't 0732 if ((typeof v) == 'string') { 0733 v = parse_number(v,undefined); 0734 if (v === undefined) return undefined; 0735 } 0736 0737 if (v != 0) { 0738 var d = new Resistor(n1,n2,v); 0739 return this.add_device(d, name); 0740 } else return this.v(n1,n2,'0',name); // zero resistance == 0V voltage source 0741 } 0742 0743 Circuit.prototype.d = function(n1,n2,area,type,name) { 0744 // try to convert string value into numeric value, barf if we can't 0745 if ((typeof area) == 'string') { 0746 area = parse_number(area,undefined); 0747 if (area === undefined) return undefined; 0748 } 0749 0750 if (area != 0) { 0751 var d = new Diode(n1,n2,area,type); 0752 return this.add_device(d, name); 0753 } // zero area diodes discarded. 0754 } 0755 0756 Circuit.prototype.c = function(n1,n2,v,name) { 0757 // try to convert string value into numeric value, barf if we can't 0758 if ((typeof v) == 'string') { 0759 v = parse_number(v,undefined); 0760 if (v === undefined) return undefined; 0761 } 0762 var d = new Capacitor(n1,n2,v); 0763 return this.add_device(d, name); 0764 } 0765 0766 Circuit.prototype.l = function(n1,n2,v,name) { 0767 // try to convert string value into numeric value, barf if we can't 0768 if ((typeof v) == 'string') { 0769 v = parse_number(v,undefined); 0770 if (v === undefined) return undefined; 0771 } 0772 var branch = this.node(undefined,T_CURRENT); 0773 var d = new Inductor(n1,n2,branch,v); 0774 return this.add_device(d, name); 0775 } 0776 0777 Circuit.prototype.v = function(n1,n2,v,name) { 0778 var branch = this.node(undefined,T_CURRENT); 0779 var d = new VSource(n1,n2,branch,v); 0780 this.voltage_sources.push(d); 0781 return this.add_device(d, name); 0782 } 0783 0784 Circuit.prototype.i = function(n1,n2,v,name) { 0785 var d = new ISource(n1,n2,v); 0786 this.current_sources.push(d); 0787 return this.add_device(d, name); 0788 } 0789 0790 Circuit.prototype.opamp = function(np,nn,no,ng,A,name) { 0791 var ratio; 0792 // try to convert string value into numeric value, barf if we can't 0793 if ((typeof A) == 'string') { 0794 ratio = parse_number(A,undefined); 0795 if (A === undefined) return undefined; 0796 } 0797 var branch = this.node(undefined,T_CURRENT); 0798 var d = new Opamp(np,nn,no,ng,branch,A,name); 0799 return this.add_device(d, name); 0800 } 0801 0802 Circuit.prototype.n = function(d,g,s, ratio, name) { 0803 // try to convert string value into numeric value, barf if we can't 0804 if ((typeof ratio) == 'string') { 0805 ratio = parse_number(ratio,undefined); 0806 if (ratio === undefined) return undefined; 0807 } 0808 var d = new Fet(d,g,s,ratio,name,'n'); 0809 return this.add_device(d, name); 0810 } 0811 0812 Circuit.prototype.p = function(d,g,s, ratio, name) { 0813 // try to convert string value into numeric value, barf if we can't 0814 if ((typeof ratio) == 'string') { 0815 ratio = parse_number(ratio,undefined); 0816 if (ratio === undefined) return undefined; 0817 } 0818 var d = new Fet(d,g,s,ratio,name,'p'); 0819 return this.add_device(d, name); 0820 } 0821 0822 /////////////////////////////////////////////////////////////////////////////// 0823 // 0824 // Support for creating conductance and capacitance matrices associated with 0825 // modified nodal analysis (unknowns are node voltages and inductor and voltage 0826 // source currents). 0827 // The linearized circuit is written as 0828 // C d/dt x = G x + rhs 0829 // x - vector of node voltages and element currents 0830 // rhs - vector of source values 0831 // C - Matrix whose values are capacitances and inductances, has many zero rows. 0832 // G - Matrix whose values are conductances and +-1's. 0833 // 0834 //////////////////////////////////////////////////////////////////////////////// 0835 0836 // add val component between two nodes to matrix M 0837 // Index of -1 refers to ground node 0838 Circuit.prototype.add_two_terminal = function(i,j,g,M) { 0839 if (i >= 0) { 0840 M[i][i] += g; 0841 if (j >= 0) { 0842 M[i][j] -= g; 0843 M[j][i] -= g; 0844 M[j][j] += g; 0845 } 0846 } else if (j >= 0) 0847 M[j][j] += g; 0848 } 0849 0850 // add val component between two nodes to matrix M 0851 // Index of -1 refers to ground node 0852 Circuit.prototype.get_two_terminal = function(i,j,x) { 0853 var xi_minus_xj = 0; 0854 if (i >= 0) xi_minus_xj = x[i]; 0855 if (j >= 0) xi_minus_xj -= x[j]; 0856 return xi_minus_xj 0857 } 0858 0859 Circuit.prototype.add_conductance_l = function(i,j,g) { 0860 this.add_two_terminal(i,j,g, this.Gl) 0861 } 0862 0863 Circuit.prototype.add_conductance = function(i,j,g) { 0864 this.add_two_terminal(i,j,g, this.G) 0865 } 0866 0867 Circuit.prototype.add_capacitance = function(i,j,c) { 0868 this.add_two_terminal(i,j,c,this.C) 0869 } 0870 0871 // add individual conductance to Gl matrix 0872 Circuit.prototype.add_to_Gl = function(i,j,g) { 0873 if (i >=0 && j >= 0) 0874 this.Gl[i][j] += g; 0875 } 0876 0877 // add individual conductance to Gl matrix 0878 Circuit.prototype.add_to_G = function(i,j,g) { 0879 if (i >=0 && j >= 0) 0880 this.G[i][j] += g; 0881 } 0882 0883 // add individual capacitance to C matrix 0884 Circuit.prototype.add_to_C = function(i,j,c) { 0885 if (i >=0 && j >= 0) 0886 this.C[i][j] += c; 0887 } 0888 0889 // add source info to rhs 0890 Circuit.prototype.add_to_rhs = function(i,v,rhs) { 0891 if (i >= 0) rhs[i] += v; 0892 } 0893 0894 0895 /////////////////////////////////////////////////////////////////////////////// 0896 // 0897 // Generic matrix support - making, copying, factoring, rank, etc 0898 // Note, Matrices are stored using nested javascript arrays. 0899 //////////////////////////////////////////////////////////////////////////////// 0900 0901 // Allocate an NxM matrix 0902 function mat_make(N,M) { 0903 var mat = new Array(N); 0904 for (var i = N - 1; i >= 0; --i) { 0905 mat[i] = new Array(M); 0906 for (var j = M - 1; j >= 0; --j) { 0907 mat[i][j] = 0.0; 0908 } 0909 } 0910 return mat; 0911 } 0912 0913 // Form b = scale*Mx 0914 function mat_v_mult(M,x,b,scale) { 0915 var n = M.length; 0916 var m = M[0].length; 0917 0918 if (n != b.length || m != x.length) 0919 throw 'Rows of M mismatched to b or cols mismatch to x.'; 0920 0921 for (var i = 0; i < n; i++) { 0922 var temp = 0; 0923 for (var j = 0; j < m; j++) temp += M[i][j]*x[j]; 0924 b[i] = scale*temp; // Recall the neg in the name 0925 } 0926 } 0927 0928 // C = scalea*A + scaleb*B, scalea, scaleb eithers numbers or arrays (row scaling) 0929 function mat_scale_add(A, B, scalea, scaleb, C) { 0930 var n = A.length; 0931 var m = A[0].length; 0932 0933 if (n > B.length || m > B[0].length) 0934 throw 'Row or columns of A to large for B'; 0935 if (n > C.length || m > C[0].length) 0936 throw 'Row or columns of A to large for C'; 0937 if ((typeof scalea == 'number') && (typeof scaleb == 'number')) 0938 for (var i = 0; i < n; i++) 0939 for (var j = 0; j < m; j++) 0940 C[i][j] = scalea*A[i][j] + scaleb*B[i][j]; 0941 else if ((typeof scaleb == 'number') && (scalea instanceof Array)) 0942 for (var i = 0; i < n; i++) 0943 for (var j = 0; j < m; j++) 0944 C[i][j] = scalea[i]*A[i][j] + scaleb*B[i][j]; 0945 else if ((typeof scaleb instanceof Array) && (scalea instanceof Array)) 0946 for (var i = 0; i < n; i++) 0947 for (var j = 0; j < m; j++) 0948 C[i][j] = scalea[i]*A[i][j] + scaleb[i]*B[i][j]; 0949 else 0950 throw 'scalea and scaleb must be scalars or Arrays'; 0951 } 0952 0953 // Returns a vector of ones and zeros, ones denote algebraic 0954 // variables (rows that can be removed without changing rank(M). 0955 Circuit.prototype.algebraic = function(M) { 0956 var Nr = M.length 0957 var Mc = mat_make(Nr, Nr); 0958 mat_copy(M,Mc); 0959 var R = mat_rank(Mc); 0960 0961 var one_if_alg = new Array(Nr); 0962 for (var row = 0; row < Nr; row++) { // psuedo gnd row small 0963 for (var col = Nr - 1; col >= 0; --col) 0964 Mc[row][col] = 0; 0965 if (mat_rank(Mc) == R) // Zeroing row left rank unchanged 0966 one_if_alg[row] = 1; 0967 else { // Zeroing row changed rank, put back 0968 for (var col = Nr - 1; col >= 0; --col) 0969 Mc[row][col] = M[row][col]; 0970 one_if_alg[row] = 0; 0971 } 0972 } 0973 return one_if_alg; 0974 } 0975 0976 // Copy A -> using the bounds of A 0977 function mat_copy(src,dest) { 0978 var n = src.length; 0979 var m = src[0].length; 0980 if (n > dest.length || m > dest[0].length) 0981 throw 'Rows or cols > rows or cols of dest'; 0982 0983 for (var i = 0; i < n; i++) 0984 for (var j = 0; j < m; j++) 0985 dest[i][j] = src[i][j]; 0986 } 0987 // Copy and transpose A -> using the bounds of A 0988 function mat_copy_transposed(src,dest) { 0989 var n = src.length; 0990 var m = src[0].length; 0991 if (n > dest[0].length || m > dest.length) 0992 throw 'Rows or cols > cols or rows of dest'; 0993 0994 for (var i = 0; i < n; i++) 0995 for (var j = 0; j < m; j++) 0996 dest[j][i] = src[i][j]; 0997 } 0998 0999 1000 // Uses GE to determine rank. 1001 function mat_rank(Mo) { 1002 var Nr = Mo.length; // Number of rows 1003 var Nc = Mo[0].length; // Number of columns 1004 var temp,i,j; 1005 // Make a copy to avoid overwriting 1006 var M = mat_make(Nr, Nc); 1007 mat_copy(Mo,M); 1008 1009 // Find matrix maximum entry 1010 var max_abs_entry = 0; 1011 for(var row = Nr-1; row >= 0; --row) { 1012 for(var col = Nr-1; col >= 0; --col) { 1013 if (Math.abs(M[row][col]) > max_abs_entry) 1014 max_abs_entry = Math.abs(M[row][col]); 1015 } 1016 } 1017 1018 // Gaussian elimination to find rank 1019 var the_rank = 0; 1020 var start_col = 0; 1021 for (var row = 0; row < Nr; row++) { 1022 // Search for first nonzero column in the remaining rows. 1023 for (var col = start_col; col < Nc; col++) { 1024 var max_v = Math.abs(M[row][col]); 1025 var max_row = row; 1026 for (var i = row + 1; i < Nr; i++) { 1027 temp = Math.abs(M[i][col]); 1028 if (temp > max_v) { max_v = temp; max_row = i; } 1029 } 1030 // if max_v non_zero, column is nonzero, eliminate in subsequent rows 1031 if (Math.abs(max_v) > eps*max_abs_entry) { 1032 start_col = col+1; 1033 the_rank += 1; 1034 // Swap rows to get max in M[row][col] 1035 temp = M[row]; 1036 M[row] = M[max_row]; 1037 M[max_row] = temp; 1038 1039 // now eliminate this column for all subsequent rows 1040 for (var i = row + 1; i < Nr; i++) { 1041 temp = M[i][col]/M[row][col]; // multiplier for current row 1042 if (temp != 0) // subtract 1043 for (var j = col; j < Nc; j++) M[i][j] -= M[row][j]*temp; 1044 } 1045 // Now move on to the next row 1046 break; 1047 } 1048 } 1049 } 1050 1051 return the_rank; 1052 } 1053 1054 // Solve Mx=b and return vector x using R^TQ^T factorization. 1055 // Multiplication by R^T implicit, should be null-space free soln. 1056 // M should have the extra column! 1057 // Almost everything is in-lined for speed, sigh. 1058 function mat_solve_rq(M, rhs) { 1059 var scale; 1060 var Nr = M.length; // Number of rows 1061 var Nc = M[0].length; // Number of columns 1062 1063 // Copy the rhs in to the last column of M if one is given. 1064 if (rhs != null) { 1065 for (var row = Nr - 1; row >= 0; --row) 1066 M[row][Nc-1] = rhs[row]; 1067 } 1068 1069 var mat_scale = 0; // Sets the scale for comparison to zero. 1070 var max_nonzero_row = Nr-1; // Assumes M nonsingular. 1071 for (var row = 0; row < Nr; row++) { 1072 // Find largest row with largest 2-norm 1073 var max_row = row; 1074 var maxsumsq = 0; 1075 for (var rowp = row; rowp < Nr; rowp++) { 1076 var Mr = M[rowp]; 1077 var sumsq = 0; 1078 for (var col = Nc-2; col >= 0; --col) // Last col=rhs 1079 sumsq += Mr[col]*Mr[col]; 1080 if ((row == rowp) || (sumsq > maxsumsq)) { 1081 max_row = rowp; 1082 maxsumsq = sumsq; 1083 } 1084 } 1085 if (max_row > row) { // Swap rows if not max row 1086 var temp = M[row]; 1087 M[row] = M[max_row]; 1088 M[max_row] = temp; 1089 } 1090 1091 // Calculate row norm, save if this is first (largest) 1092 var row_norm = Math.sqrt(maxsumsq); 1093 if (row == 0) mat_scale = row_norm; 1094 1095 // Check for all zero rows 1096 if (row_norm > mat_scale*eps) 1097 scale = 1.0/row_norm; 1098 else { 1099 max_nonzero_row = row - 1; // Rest will be nullspace of M 1100 break; 1101 } 1102 1103 // Nonzero row, eliminate from rows below 1104 var Mr = M[row]; 1105 for (var col = Nc-1; col >= 0; --col) // Scale rhs also 1106 Mr[col] *= scale; 1107 for (var rowp = row + 1; rowp < Nr; rowp++) { // Update. 1108 var Mrp = M[rowp]; 1109 var inner = 0; 1110 for (var col = Nc-2; col >= 0; --col) // Project 1111 inner += Mr[col]*Mrp[col]; 1112 for (var col = Nc-1; col >= 0; --col) // Ortho (rhs also) 1113 Mrp[col] -= inner *Mr[col]; 1114 } 1115 } 1116 1117 // Last Column of M has inv(R^T)*rhs. Scale rows of Q to get x. 1118 var x = new Array(Nc-1); 1119 for (var col = Nc-2; col >= 0; --col) 1120 x[col] = 0; 1121 for (var row = max_nonzero_row; row >= 0; --row) { 1122 Mr = M[row]; 1123 for (var col = Nc-2; col >= 0; --col) { 1124 x[col] += Mr[col]*Mr[Nc-1]; 1125 } 1126 } 1127 1128 return x; 1129 } 1130 1131 // solve Mx=b and return vector x given augmented matrix M = [A | b] 1132 // Uses Gaussian elimination with partial pivoting 1133 function mat_solve(M,rhs) { 1134 var N = M.length; // augmented matrix M has N rows, N+1 columns 1135 var temp,i,j; 1136 1137 // Copy the rhs in to the last column of M if one is given. 1138 if (rhs != null) { 1139 for (var row = 0; row < N ; row++) 1140 M[row][N] = rhs[row]; 1141 } 1142 1143 // gaussian elimination 1144 for (var col = 0; col < N ; col++) { 1145 // find pivot: largest abs(v) in this column of remaining rows 1146 var max_v = Math.abs(M[col][col]); 1147 var max_col = col; 1148 for (i = col + 1; i < N; i++) { 1149 temp = Math.abs(M[i][col]); 1150 if (temp > max_v) { max_v = temp; max_col = i; } 1151 } 1152 1153 // if no value found, generate a small conductance to gnd 1154 // otherwise swap current row with pivot row 1155 if (max_v == 0) M[col][col] = eps; 1156 else { 1157 temp = M[col]; 1158 M[col] = M[max_col]; 1159 M[max_col] = temp; 1160 } 1161 1162 // now eliminate this column for all subsequent rows 1163 for (i = col + 1; i < N; i++) { 1164 temp = M[i][col]/M[col][col]; // multiplier we'll use for current row 1165 if (temp != 0) 1166 // subtract current row from row we're working on 1167 // remember to process b too! 1168 for (j = col; j <= N; j++) M[i][j] -= M[col][j]*temp; 1169 } 1170 } 1171 1172 // matrix is now upper triangular, so solve for elements of x starting 1173 // with the last row 1174 var x = new Array(N); 1175 for (i = N-1; i >= 0; --i) { 1176 temp = M[i][N]; // grab b[i] from augmented matrix as RHS 1177 // subtract LHS term from RHS using known x values 1178 for (j = N-1; j > i; --j) temp -= M[i][j]*x[j]; 1179 // now compute new x value 1180 x[i] = temp/M[i][i]; 1181 } 1182 1183 return x; 1184 } 1185 1186 // test solution code, expect x = [2,3,-1] 1187 //M = [[2,1,-1,8],[-3,-1,2,-11],[-2,1,2,-3]]; 1188 //x = mat_solve(M); 1189 //y = 1; // so we have place to set a breakpoint :) 1190 1191 /////////////////////////////////////////////////////////////////////////////// 1192 // 1193 // Device base class 1194 // 1195 //////////////////////////////////////////////////////////////////////////////// 1196 1197 function Device() { 1198 } 1199 1200 // complete initial set up of device 1201 Device.prototype.finalize = function() { 1202 } 1203 1204 // Load the linear elements in to Gl and C 1205 Device.prototype.load_linear = function(ckt) { 1206 } 1207 1208 // load linear system equations for dc analysis 1209 // (inductors shorted and capacitors opened) 1210 Device.prototype.load_dc = function(ckt,soln,rhs) { 1211 } 1212 1213 // load linear system equations for tran analysis 1214 Device.prototype.load_tran = function(ckt,soln) { 1215 } 1216 1217 // load linear system equations for ac analysis: 1218 // current sources open, voltage sources shorted 1219 // linear models at operating point for everyone else 1220 Device.prototype.load_ac = function(ckt,rhs) { 1221 } 1222 1223 // return time of next breakpoint for the device 1224 Device.prototype.breakpoint = function(time) { 1225 return undefined; 1226 } 1227 1228 /////////////////////////////////////////////////////////////////////////////// 1229 // 1230 // Parse numbers in engineering notation 1231 // 1232 /////////////////////////////////////////////////////////////////////////////// 1233 1234 // convert first character of argument into an integer 1235 function ord(ch) { 1236 return ch.charCodeAt(0); 1237 } 1238 1239 // convert string argument to a number, accepting usual notations 1240 // (hex, octal, binary, decimal, floating point) plus engineering 1241 // scale factors (eg, 1k = 1000.0 = 1e3). 1242 // return default if argument couldn't be interpreted as a number 1243 function parse_number(s,default_v) { 1244 var slen = s.length; 1245 var multiplier = 1; 1246 var result = 0; 1247 var index = 0; 1248 1249 // skip leading whitespace 1250 while (index < slen && s.charAt(index) <= ' ') index += 1; 1251 if (index == slen) return default_v; 1252 1253 // check for leading sign 1254 if (s.charAt(index) == '-') { 1255 multiplier = -1; 1256 index += 1; 1257 } else if (s.charAt(index) == '+') 1258 index += 1; 1259 var start = index; // remember where digits start 1260 1261 // if leading digit is 0, check for hex, octal or binary notation 1262 if (index >= slen) return default_v; 1263 else if (s.charAt(index) == '0') { 1264 index += 1; 1265 if (index >= slen) return 0; 1266 if (s.charAt(index) == 'x' || s.charAt(index) == 'X') { // hex 1267 while (true) { 1268 index += 1; 1269 if (index >= slen) break; 1270 if (s.charAt(index) >= '0' && s.charAt(index) <= '9') 1271 result = result*16 + ord(s.charAt(index)) - ord('0'); 1272 else if (s.charAt(index) >= 'A' && s.charAt(index) <= 'F') 1273 result = result*16 + ord(s.charAt(index)) - ord('A') + 10; 1274 else if (s.charAt(index) >= 'a' && s.charAt(index) <= 'f') 1275 result = result*16 + ord(s.charAt(index)) - ord('a') + 10; 1276 else break; 1277 } 1278 return result*multiplier; 1279 } else if (s.charAt(index) == 'b' || s.charAt(index) == 'B') { // binary 1280 while (true) { 1281 index += 1; 1282 if (index >= slen) break; 1283 if (s.charAt(index) >= '0' && s.charAt(index) <= '1') 1284 result = result*2 + ord(s.charAt(index)) - ord('0'); 1285 else break; 1286 } 1287 return result*multiplier; 1288 } else if (s.charAt(index) != '.') { // octal 1289 while (true) { 1290 if (s.charAt(index) >= '0' && s.charAt(index) <= '7') 1291 result = result*8 + ord(s.charAt(index)) - ord('0'); 1292 else break; 1293 index += 1; 1294 if (index >= slen) break; 1295 } 1296 return result*multiplier; 1297 } 1298 } 1299 // read decimal integer or floating-point number 1300 while (true) { 1301 if (s.charAt(index) >= '0' && s.charAt(index) <= '9') 1302 result = result*10 + ord(s.charAt(index)) - ord('0'); 1303 else break; 1304 index += 1; 1305 if (index >= slen) break; 1306 } 1307 1308 // fractional part? 1309 if (index < slen && s.charAt(index) == '.') { 1310 while (true) { 1311 index += 1; 1312 if (index >= slen) break; 1313 if (s.charAt(index) >= '0' && s.charAt(index) <= '9') { 1314 result = result*10 + ord(s.charAt(index)) - ord('0'); 1315 multiplier *= 0.1; 1316 } else break; 1317 } 1318 } 1319 1320 // if we haven't seen any digits yet, don't check 1321 // for exponents or scale factors 1322 if (index == start) return default_v; 1323 1324 // type of multiplier determines type of result: 1325 // multiplier is a float if we've seen digits past 1326 // a decimal point, otherwise it's an int or long. 1327 // Up to this point result is an int or long. 1328 result *= multiplier; 1329 1330 // now check for exponent or engineering scale factor. If there 1331 // is one, result will be a float. 1332 if (index < slen) { 1333 var scale = s.charAt(index); 1334 index += 1; 1335 if (scale == 'e' || scale == 'E') { 1336 var exponent = 0; 1337 multiplier = 10.0; 1338 if (index < slen) { 1339 if (s.charAt(index) == '+') index += 1; 1340 else if (s.charAt(index) == '-') { 1341 index += 1; 1342 multiplier = 0.1; 1343 } 1344 } 1345 while (index < slen) { 1346 if (s.charAt(index) >= '0' && s.charAt(index) <= '9') { 1347 exponent = exponent*10 + ord(s.charAt(index)) - ord('0'); 1348 index += 1; 1349 } else break; 1350 } 1351 while (exponent > 0) { 1352 exponent -= 1; 1353 result *= multiplier; 1354 } 1355 } else if (scale == 't' || scale == 'T') result *= 1e12; 1356 else if (scale == 'g' || scale == 'G') result *= 1e9; 1357 else if (scale == 'M') result *= 1e6; 1358 else if (scale == 'k' || scale == 'K') result *= 1e3; 1359 else if (scale == 'm') result *= 1e-3; 1360 else if (scale == 'u' || scale == 'U') result *= 1e-6; 1361 else if (scale == 'n' || scale == 'N') result *= 1e-9; 1362 else if (scale == 'p' || scale == 'P') result *= 1e-12; 1363 else if (scale == 'f' || scale == 'F') result *= 1e-15; 1364 } 1365 // ignore any remaining chars, eg, 1kohms returns 1000 1366 return result; 1367 } 1368 1369 Circuit.prototype.parse_number = parse_number; // make it easy to call from outside 1370 1371 /////////////////////////////////////////////////////////////////////////////// 1372 // 1373 // Sources 1374 // 1375 /////////////////////////////////////////////////////////////////////////////// 1376 1377 // argument is a string describing the source's value (see comments for details) 1378 // source types: dc,step,square,triangle,sin,pulse,pwl,pwl_repeating 1379 1380 // returns an object with the following attributes: 1381 // fun -- name of source function 1382 // args -- list of argument values 1383 // value(t) -- compute source value at time t 1384 // inflection_point(t) -- compute time after t when a time point is needed 1385 // dc -- value at time 0 1386 // period -- repeat period for periodic sources (0 if not periodic) 1387 1388 function parse_source(v) { 1389 // generic parser: parse v as either <value> or <fun>(<value>,...) 1390 var src = {}; 1391 src.period = 0; // Default not periodic 1392 src.value = function(t) { return 0; } // overridden below 1393 src.inflection_point = function(t) { return undefined; }; // may be overridden below 1394 1395 // see if there's a "(" in the description 1396 var index = v.indexOf('('); 1397 var ch; 1398 if (index >= 0) { 1399 src.fun = v.slice(0,index); // function name is before the "(" 1400 src.args = []; // we'll push argument values onto this list 1401 var end = v.indexOf(')',index); 1402 if (end == -1) end = v.length; 1403 1404 index += 1; // start parsing right after "(" 1405 while (index < end) { 1406 // figure out where next argument value starts 1407 ch = v.charAt(index); 1408 if (ch <= ' ') { index++; continue; } 1409 // and where it ends 1410 var arg_end = v.indexOf(',',index); 1411 if (arg_end == -1) arg_end = end; 1412 // parse and save result in our list of arg values 1413 src.args.push(parse_number(v.slice(index,arg_end),undefined)); 1414 index = arg_end + 1; 1415 } 1416 } else { 1417 src.fun = 'dc'; 1418 src.args = [parse_number(v,0)]; 1419 } 1420 1421 // post-processing for constant sources 1422 // dc(v) 1423 if (src.fun == 'dc') { 1424 var v = arg_value(src.args,0,0); 1425 src.args = [v]; 1426 src.value = function(t) { return v; } // closure 1427 } 1428 1429 // post-processing for impulse sources 1430 // impulse(height,width) 1431 else if (src.fun == 'impulse') { 1432 var h = arg_value(src.args,0,1); // default height: 1 1433 var w = Math.abs(arg_value(src.args,2,1e-9)); // default width: 1ns 1434 src.args = [h,w]; // remember any defaulted values 1435 pwl_source(src,[0,0,w/2,h,w,0],false); 1436 } 1437 1438 // post-processing for step sources 1439 // step(v_init,v_plateau,t_delay,t_rise) 1440 else if (src.fun == 'step') { 1441 var v1 = arg_value(src.args,0,0); // default init value: 0V 1442 var v2 = arg_value(src.args,1,1); // default plateau value: 1V 1443 var td = Math.max(0,arg_value(src.args,2,0)); // time step starts 1444 var tr = Math.abs(arg_value(src.args,3,1e-9)); // default rise time: 1ns 1445 src.args = [v1,v2,td,tr]; // remember any defaulted values 1446 pwl_source(src,[td,v1,td+tr,v2],false); 1447 } 1448 1449 // post-processing for square wave 1450 // square(v_init,v_plateau,freq,duty_cycle) 1451 else if (src.fun == 'square') { 1452 var v1 = arg_value(src.args,0,0); // default init value: 0V 1453 var v2 = arg_value(src.args,1,1); // default plateau value: 1V 1454 var freq = Math.abs(arg_value(src.args,2,1)); // default frequency: 1Hz 1455 var duty_cycle = Math.min(100,Math.abs(arg_value(src.args,3,50))); // default duty cycle: 0.5 1456 src.args = [v1,v2,freq,duty_cycle]; // remember any defaulted values 1457 1458 var per = freq == 0 ? Infinity : 1/freq; 1459 var t_change = 0.01 * per; // rise and fall time 1460 var t_pw = .01 * duty_cycle * 0.98 * per; // fraction of cycle minus rise and fall time 1461 pwl_source(src,[0,v1,t_change,v2,t_change+t_pw, 1462 v2,t_change+t_pw+t_change,v1,per,v1],true); 1463 } 1464 1465 // post-processing for triangle 1466 // triangle(v_init,v_plateua,t_period) 1467 else if (src.fun == 'triangle') { 1468 var v1 = arg_value(src.args,0,0); // default init value: 0V 1469 var v2 = arg_value(src.args,1,1); // default plateau value: 1V 1470 var freq = Math.abs(arg_value(src.args,2,1)); // default frequency: 1s 1471 src.args = [v1,v2,freq]; // remember any defaulted values 1472 1473 var per = freq == 0 ? Infinity : 1/freq; 1474 pwl_source(src,[0,v1,per/2,v2,per,v1],true); 1475 } 1476 1477 // post-processing for pwl and pwlr sources 1478 // pwl[r](t1,v1,t2,v2,...) 1479 else if (src.fun == 'pwl' || src.fun == 'pwl_repeating') { 1480 pwl_source(src,src.args,src.fun == 'pwl_repeating'); 1481 } 1482 1483 // post-processing for pulsed sources 1484 // pulse(v_init,v_plateau,t_delay,t_rise,t_fall,t_width,t_period) 1485 else if (src.fun == 'pulse') { 1486 var v1 = arg_value(src.args,0,0); // default init value: 0V 1487 var v2 = arg_value(src.args,1,1); // default plateau value: 1V 1488 var td = Math.max(0,arg_value(src.args,2,0)); // time pulse starts 1489 var tr = Math.abs(arg_value(src.args,3,1e-9)); // default rise time: 1ns 1490 var tf = Math.abs(arg_value(src.args,4,1e-9)); // default rise time: 1ns 1491 var pw = Math.abs(arg_value(src.args,5,1e9)); // default pulse width: "infinite" 1492 var per = Math.abs(arg_value(src.args,6,1e9)); // default period: "infinite" 1493 src.args = [v1,v2,td,tr,tf,pw,per]; 1494 1495 var t1 = td; // time when v1 -> v2 transition starts 1496 var t2 = t1 + tr; // time when v1 -> v2 transition ends 1497 var t3 = t2 + pw; // time when v2 -> v1 transition starts 1498 var t4 = t3 + tf; // time when v2 -> v1 transition ends 1499 1500 pwl_source(src,[t1,v1, t2,v2, t3,v2, t4,v1, per,v1],true); 1501 } 1502 1503 // post-processing for sinusoidal sources 1504 // sin(v_offset,v_amplitude,freq_hz,t_delay,phase_offset_degrees) 1505 else if (src.fun == 'sin') { 1506 var voffset = arg_value(src.args,0,0); // default offset voltage: 0V 1507 var va = arg_value(src.args,1,1); // default amplitude: -1V to 1V 1508 var freq = Math.abs(arg_value(src.args,2,1)); // default frequency: 1Hz 1509 src.period = 1.0/freq; 1510 1511 var td = Math.max(0,arg_value(src.args,3,0)); // default time delay: 0sec 1512 var phase = arg_value(src.args,4,0); // default phase offset: 0 degrees 1513 src.args = [voffset,va,freq,td,phase]; 1514 1515 phase /= 360.0; 1516 1517 // return value of source at time t 1518 src.value = function(t) { // closure 1519 if (t < td) return voffset + va*Math.sin(2*Math.PI*phase); 1520 else return voffset + va*Math.sin(2*Math.PI*(freq*(t - td) + phase)); 1521 } 1522 1523 // return time of next inflection point after time t 1524 src.inflection_point = function(t) { // closure 1525 if (t < td) return td; 1526 else return undefined; 1527 } 1528 } 1529 1530 // object has all the necessary info to compute the source value and inflection points 1531 src.dc = src.value(0); // DC value is value at time 0 1532 return src; 1533 } 1534 1535 function pwl_source(src,tv_pairs,repeat) { 1536 var nvals = tv_pairs.length; 1537 if (repeat) 1538 src.period = tv_pairs[nvals-2]; // Repeat period of source 1539 if (nvals % 2 == 1) npts -= 1; // make sure it's even! 1540 1541 if (nvals <= 2) { 1542 // handle degenerate case 1543 src.value = function(t) { return nvals == 2 ? tv_pairs[1] : 0; } 1544 src.inflection_point = function(t) { return undefined; } 1545 } else { 1546 src.value = function(t) { // closure 1547 if (repeat) 1548 // make time periodic if values are to be repeated 1549 t = Math.fmod(t,tv_pairs[nvals-2]); 1550 var last_t = tv_pairs[0]; 1551 var last_v = tv_pairs[1]; 1552 if (t > last_t) { 1553 var next_t,next_v; 1554 for (var i = 2; i < nvals; i += 2) { 1555 next_t = tv_pairs[i]; 1556 next_v = tv_pairs[i+1]; 1557 if (next_t > last_t) // defend against bogus tv pairs 1558 if (t < next_t) 1559 return last_v + (next_v - last_v)*(t - last_t)/(next_t - last_t); 1560 last_t = next_t; 1561 last_v = next_v; 1562 } 1563 } 1564 return last_v; 1565 } 1566 src.inflection_point = function(t) { // closure 1567 if (repeat) 1568 // make time periodic if values are to be repeated 1569 t = Math.fmod(t,tv_pairs[nvals-2]); 1570 for (var i = 0; i < nvals; i += 2) { 1571 var next_t = tv_pairs[i]; 1572 if (t < next_t) return next_t; 1573 } 1574 return undefined; 1575 } 1576 } 1577 } 1578 1579 // helper function: return args[index] if present, else default_v 1580 function arg_value(args,index,default_v) { 1581 if (index < args.length) { 1582 var result = args[index]; 1583 if (result === undefined) result = default_v; 1584 return result; 1585 } else return default_v; 1586 } 1587 1588 // we need fmod in the Math library! 1589 Math.fmod = function(numerator,denominator) { 1590 var quotient = Math.floor(numerator/denominator); 1591 return numerator - quotient*denominator; 1592 } 1593 1594 /////////////////////////////////////////////////////////////////////////////// 1595 // 1596 // Sources 1597 // 1598 /////////////////////////////////////////////////////////////////////////////// 1599 1600 function VSource(npos,nneg,branch,v) { 1601 Device.call(this); 1602 this.src = parse_source(v); 1603 this.npos = npos; 1604 this.nneg = nneg; 1605 this.branch = branch; 1606 } 1607 VSource.prototype = new Device(); 1608 VSource.prototype.constructor = VSource; 1609 1610 // load linear part for source evaluation 1611 VSource.prototype.load_linear = function(ckt) { 1612 // MNA stamp for independent voltage source 1613 ckt.add_to_Gl(this.branch,this.npos,1.0); 1614 ckt.add_to_Gl(this.branch,this.nneg,-1.0); 1615 ckt.add_to_Gl(this.npos,this.branch,1.0); 1616 ckt.add_to_Gl(this.nneg,this.branch,-1.0); 1617 } 1618 1619 // Source voltage added to b. 1620 VSource.prototype.load_dc = function(ckt,soln,rhs) { 1621 ckt.add_to_rhs(this.branch,this.src.dc,rhs); 1622 } 1623 1624 // Load time-dependent value for voltage source for tran 1625 VSource.prototype.load_tran = function(ckt,soln,rhs,time) { 1626 ckt.add_to_rhs(this.branch,this.src.value(time),rhs); 1627 } 1628 1629 // return time of next breakpoint for the device 1630 VSource.prototype.breakpoint = function(time) { 1631 return this.src.inflection_point(time); 1632 } 1633 1634 // small signal model ac value 1635 VSource.prototype.load_ac = function(ckt,rhs) { 1636 ckt.add_to_rhs(this.branch,1.0,rhs); 1637 } 1638 1639 function ISource(npos,nneg,v) { 1640 Device.call(this); 1641 this.src = parse_source(v); 1642 this.npos = npos; 1643 this.nneg = nneg; 1644 } 1645 ISource.prototype = new Device(); 1646 ISource.prototype.constructor = ISource; 1647 1648 ISource.prototype.load_linear = function(ckt) { 1649 // Current source is open when off, no linear contribution 1650 } 1651 1652 // load linear system equations for dc analysis 1653 ISource.prototype.load_dc = function(ckt,soln,rhs) { 1654 var is = this.src.dc; 1655 1656 // MNA stamp for independent current source 1657 ckt.add_to_rhs(this.npos,-is,rhs); // current flow into npos 1658 ckt.add_to_rhs(this.nneg,is,rhs); // and out of nneg 1659 } 1660 1661 // load linear system equations for tran analysis (just like DC) 1662 ISource.prototype.load_tran = function(ckt,soln,rhs,time) { 1663 var is = this.src.value(time); 1664 1665 // MNA stamp for independent current source 1666 ckt.add_to_rhs(this.npos,-is,rhs); // current flow into npos 1667 ckt.add_to_rhs(this.nneg,is,rhs); // and out of nneg 1668 } 1669 1670 // return time of next breakpoint for the device 1671 ISource.prototype.breakpoint = function(time) { 1672 return this.src.inflection_point(time); 1673 } 1674 1675 // small signal model: open circuit 1676 ISource.prototype.load_ac = function(ckt,rhs) { 1677 // MNA stamp for independent current source 1678 ckt.add_to_rhs(this.npos,-1.0,rhs); // current flow into npos 1679 ckt.add_to_rhs(this.nneg,1.0,rhs); // and out of nneg 1680 } 1681 1682 /////////////////////////////////////////////////////////////////////////////// 1683 // 1684 // Resistor 1685 // 1686 /////////////////////////////////////////////////////////////////////////////// 1687 1688 function Resistor(n1,n2,v) { 1689 Device.call(this); 1690 this.n1 = n1; 1691 this.n2 = n2; 1692 this.g = 1.0/v; 1693 } 1694 Resistor.prototype = new Device(); 1695 Resistor.prototype.constructor = Resistor; 1696 1697 Resistor.prototype.load_linear = function(ckt) { 1698 // MNA stamp for admittance g 1699 ckt.add_conductance_l(this.n1,this.n2,this.g); 1700 } 1701 1702 Resistor.prototype.load_dc = function(ckt) { 1703 // Nothing to see here, move along. 1704 } 1705 1706 Resistor.prototype.load_tran = function(ckt,soln) { 1707 } 1708 1709 Resistor.prototype.load_ac = function(ckt) { 1710 } 1711 1712 /////////////////////////////////////////////////////////////////////////////// 1713 // 1714 // Diode 1715 // 1716 /////////////////////////////////////////////////////////////////////////////// 1717 1718 function Diode(n1,n2,v,type) { 1719 Device.call(this); 1720 this.anode = n1; 1721 this.cathode = n2; 1722 this.area = v; 1723 this.type = type; // 'normal' or 'ideal' 1724 this.is = 1.0e-14; 1725 this.ais = this.area * this.is; 1726 this.vt = (type == 'normal') ? 25.8e-3 : 0.1e-3; // 26mv or .1mv 1727 this.exp_arg_max = 50; // less than single precision max. 1728 this.exp_max = Math.exp(this.exp_arg_max); 1729 } 1730 Diode.prototype = new Device(); 1731 Diode.prototype.constructor = Diode; 1732 1733 Diode.prototype.load_linear = function(ckt) { 1734 // Diode is not linear, has no linear piece. 1735 } 1736 1737 Diode.prototype.load_dc = function(ckt,soln,rhs) { 1738 var vd = ckt.get_two_terminal(this.anode, this.cathode, soln); 1739 var exp_arg = vd / this.vt; 1740 var temp1, temp2; 1741 // Estimate exponential with a quadratic if arg too big. 1742 var abs_exp_arg = Math.abs(exp_arg); 1743 var d_arg = abs_exp_arg - this.exp_arg_max; 1744 if (d_arg > 0) { 1745 var quad = 1 + d_arg + 0.5*d_arg*d_arg; 1746 temp1 = this.exp_max * quad; 1747 temp2 = this.exp_max * (1 + d_arg); 1748 } else { 1749 temp1 = Math.exp(abs_exp_arg); 1750 temp2 = temp1; 1751 } 1752 if (exp_arg < 0) { // Use exp(-x) = 1.0/exp(x) 1753 temp1 = 1.0/temp1; 1754 temp2 = (temp1*temp2)*temp1; 1755 } 1756 var id = this.ais * (temp1 - 1); 1757 var gd = this.ais * (temp2 / this.vt); 1758 1759 // MNA stamp for independent current source 1760 ckt.add_to_rhs(this.anode,-id,rhs); // current flows into anode 1761 ckt.add_to_rhs(this.cathode,id,rhs); // and out of cathode 1762 ckt.add_conductance(this.anode,this.cathode,gd); 1763 } 1764 1765 Diode.prototype.load_tran = function(ckt,soln,rhs,time) { 1766 this.load_dc(ckt,soln,rhs); 1767 } 1768 1769 Diode.prototype.load_ac = function(ckt) { 1770 } 1771 1772 1773 /////////////////////////////////////////////////////////////////////////////// 1774 // 1775 // Capacitor 1776 // 1777 /////////////////////////////////////////////////////////////////////////////// 1778 1779 function Capacitor(n1,n2,v) { 1780 Device.call(this); 1781 this.n1 = n1; 1782 this.n2 = n2; 1783 this.value = v; 1784 } 1785 Capacitor.prototype = new Device(); 1786 Capacitor.prototype.constructor = Capacitor; 1787 1788 Capacitor.prototype.load_linear = function(ckt) { 1789 // MNA stamp for capacitance matrix 1790 ckt.add_capacitance(this.n1,this.n2,this.value); 1791 } 1792 1793 Capacitor.prototype.load_dc = function(ckt,soln,rhs) { 1794 } 1795 1796 Capacitor.prototype.load_ac = function(ckt) { 1797 } 1798 1799 Capacitor.prototype.load_tran = function(ckt) { 1800 } 1801 1802 /////////////////////////////////////////////////////////////////////////////// 1803 // 1804 // Inductor 1805 // 1806 /////////////////////////////////////////////////////////////////////////////// 1807 1808 function Inductor(n1,n2,branch,v) { 1809 Device.call(this); 1810 this.n1 = n1; 1811 this.n2 = n2; 1812 this.branch = branch; 1813 this.value = v; 1814 } 1815 Inductor.prototype = new Device(); 1816 Inductor.prototype.constructor = Inductor; 1817 1818 Inductor.prototype.load_linear = function(ckt) { 1819 // MNA stamp for inductor linear part 1820 // L on diag of C because L di/dt = v(n1) - v(n2) 1821 ckt.add_to_Gl(this.n1,this.branch,1); 1822 ckt.add_to_Gl(this.n2,this.branch,-1); 1823 ckt.add_to_Gl(this.branch,this.n1,-1); 1824 ckt.add_to_Gl(this.branch,this.n2,1); 1825 ckt.add_to_C(this.branch,this.branch,this.value) 1826 } 1827 1828 Inductor.prototype.load_dc = function(ckt,soln,rhs) { 1829 // Inductor is a short at dc, so is linear. 1830 } 1831 1832 Inductor.prototype.load_ac = function(ckt) { 1833 } 1834 1835 Inductor.prototype.load_tran = function(ckt) { 1836 } 1837 1838 1839 /////////////////////////////////////////////////////////////////////////////// 1840 // 1841 // Simple Voltage-Controlled Voltage Source Op Amp model 1842 // 1843 /////////////////////////////////////////////////////////////////////////////// 1844 1845 function Opamp(np,nn,no,ng,branch,A,name) { 1846 Device.call(this); 1847 this.np = np; 1848 this.nn = nn; 1849 this.no = no; 1850 this.ng = ng; 1851 this.branch = branch; 1852 this.gain = A; 1853 this.name = name; 1854 } 1855 1856 Opamp.prototype = new Device(); 1857 Opamp.prototype.constructor = Opamp; 1858 Opamp.prototype.load_linear = function(ckt) { 1859 // MNA stamp for VCVS: 1/A(v(no) - v(ng)) - (v(np)-v(nn))) = 0. 1860 var invA = 1.0/this.gain; 1861 ckt.add_to_Gl(this.no,this.branch,1); 1862 ckt.add_to_Gl(this.ng,this.branch,-1); 1863 ckt.add_to_Gl(this.branch,this.no,invA); 1864 ckt.add_to_Gl(this.branch,this.ng,-invA); 1865 ckt.add_to_Gl(this.branch,this.np,-1); 1866 ckt.add_to_Gl(this.branch,this.nn,1); 1867 } 1868 1869 Opamp.prototype.load_dc = function(ckt,soln,rhs) { 1870 // Op-amp is linear. 1871 } 1872 1873 Opamp.prototype.load_ac = function(ckt) { 1874 } 1875 1876 Opamp.prototype.load_tran = function(ckt) { 1877 } 1878 1879 1880 /////////////////////////////////////////////////////////////////////////////// 1881 // 1882 // Simplified MOS FET with no bulk connection and no body effect. 1883 // 1884 /////////////////////////////////////////////////////////////////////////////// 1885 1886 function Fet(d,g,s,ratio,name,type) { 1887 Device.call(this); 1888 this.d = d; 1889 this.g = g; 1890 this.s = s; 1891 this.name = name; 1892 this.ratio = ratio; 1893 if (type != 'n' && type != 'p') 1894 { throw 'fet type is not n or p'; 1895 } 1896 this.type_sign = (type == 'n') ? 1 : -1; 1897 this.vt = 0.5; 1898 this.kp = 20e-6; 1899 this.beta = this.kp * this.ratio; 1900 this.lambda = 0.05; 1901 } 1902 Fet.prototype = new Device(); 1903 Fet.prototype.constructor = Fet; 1904 1905 Fet.prototype.load_linear = function(ckt) { 1906 // FET's are nonlinear, just like javascript progammers 1907 } 1908 1909 Fet.prototype.load_dc = function(ckt,soln,rhs) { 1910 var vds = this.type_sign * ckt.get_two_terminal(this.d, this.s, soln); 1911 if (vds < 0) { // Drain and source have swapped roles 1912 var temp = this.d; 1913 this.d = this.s; 1914 this.s = temp; 1915 vds = this.type_sign * ckt.get_two_terminal(this.d, this.s, soln); 1916 } 1917 var vgs = this.type_sign * ckt.get_two_terminal(this.g, this.s, soln); 1918 var vgst = vgs - this.vt; 1919 var gmgs,ids,gds; 1920 if (vgst > 0.0 ) { // vgst < 0, transistor off, no subthreshold here. 1921 if (vgst < vds) { /* Saturation. */ 1922 gmgs = this.beta * (1 + (this.lambda * vds)) * vgst; 1923 ids = this.type_sign * 0.5 * gmgs * vgst; 1924 gds = 0.5 * this.beta * vgst * vgst * this.lambda; 1925 } else { /* Linear region */ 1926 gmgs = this.beta * (1 + this.lambda * vds); 1927 ids = this.type_sign * gmgs * vds * (vgst - 0.50 * vds); 1928 gds = gmgs * (vgst - vds) + this.beta * this.lambda * vds * (vgst - 0.5 * vds); 1929 gmgs *= vds; 1930 } 1931 ckt.add_to_rhs(this.d,-ids,rhs); // current flows into the drain 1932 ckt.add_to_rhs(this.s, ids,rhs); // and out the source 1933 ckt.add_conductance(this.d,this.s,gds); 1934 ckt.add_to_G(this.s,this.s, gmgs); 1935 ckt.add_to_G(this.d,this.s,-gmgs); 1936 ckt.add_to_G(this.d,this.g, gmgs); 1937 ckt.add_to_G(this.s,this.g,-gmgs); 1938 } 1939 } 1940 1941 Fet.prototype.load_tran = function(ckt,soln,rhs) { 1942 this.load_dc(ckt,soln,rhs); 1943 } 1944 1945 Fet.prototype.load_ac = function(ckt) { 1946 } 1947 1948 1949 /////////////////////////////////////////////////////////////////////////////// 1950 // 1951 // Module definition 1952 // 1953 /////////////////////////////////////////////////////////////////////////////// 1954 var module = { 1955 'Circuit': Circuit, 1956 'parse_number': parse_number, 1957 'parse_source': parse_source 1958 } 1959 return module; 1960 }());