Qucs-core
0.0.19
|
00001 /* 00002 * diode.cpp - diode class implementation 00003 * 00004 * Copyright (C) 2004, 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include "component.h" 00030 #include "device.h" 00031 #include "devstates.h" 00032 #include "diode.h" 00033 00034 #define NODE_C 0 /* cathode node */ 00035 #define NODE_A 1 /* anode node */ 00036 00037 #define StateVars 1 // state variables 00038 00039 // state variable indices 00040 #define _UdPrev 0 00041 00042 // state variable shortcuts 00043 #define UdPrev deviceVar (_UdPrev) 00044 00045 using namespace qucs; 00046 using namespace qucs::device; 00047 00048 // Constructor for the diode. 00049 diode::diode () : circuit (2) { 00050 rs = NULL; 00051 type = CIR_DIODE; 00052 } 00053 00054 // Callback for S-parameter analysis. 00055 void diode::calcSP (nr_double_t frequency) { 00056 nr_double_t gd = getOperatingPoint ("gd"); 00057 nr_double_t Cd = getOperatingPoint ("Cd"); 00058 nr_complex_t y = 2 * z0 * nr_complex_t (gd, Cd * 2.0 * pi * frequency); 00059 setS (NODE_C, NODE_C, 1.0 / (1.0 + y)); 00060 setS (NODE_A, NODE_A, 1.0 / (1.0 + y)); 00061 setS (NODE_C, NODE_A, y / (1.0 + y)); 00062 setS (NODE_A, NODE_C, y / (1.0 + y)); 00063 } 00064 00065 // Callback for S-parameter noise analysis. 00066 void diode::calcNoiseSP (nr_double_t frequency) { 00067 #if MICHAEL /* shot noise only */ 00068 nr_double_t Id = getOperatingPoint ("Id"); 00069 nr_double_t Is = getPropertyDouble ("Is") + getPropertyDouble ("Isr"); 00070 00071 // adjust shot noise current if necessary 00072 if (Id < -Is) Id = -Is; 00073 00074 nr_double_t gd = getOperatingPoint ("gd"); 00075 nr_double_t Cd = getOperatingPoint ("Cd"); 00076 00077 nr_complex_t y = rect (gd, Cd * 2.0 * pi * frequency); 00078 nr_complex_t f = 2 * z0 * (Id + 2 * Is) / norm (2 * z0 * y + 1) * QoverkB / T0; 00079 setN (NODE_C, NODE_C, +f); setN (NODE_A, NODE_A, +f); 00080 setN (NODE_C, NODE_A, -f); setN (NODE_A, NODE_C, -f); 00081 #else 00082 setMatrixN (cytocs (calcMatrixCy (frequency) * z0, getMatrixS ())); 00083 #endif 00084 } 00085 00086 // Computes noise correlation matrix Cy. 00087 matrix diode::calcMatrixCy (nr_double_t frequency) { 00088 // fetch computed operating points 00089 nr_double_t Id = getOperatingPoint ("Id"); 00090 nr_double_t Is = getPropertyDouble ("Is") + getPropertyDouble ("Isr"); 00091 00092 // adjust shot noise current if necessary 00093 if (Id < -Is) Id = -Is; 00094 00095 nr_double_t Kf = getPropertyDouble ("Kf"); 00096 nr_double_t Af = getPropertyDouble ("Af"); 00097 nr_double_t Ffe = getPropertyDouble ("Ffe"); 00098 00099 // build noise current correlation matrix 00100 matrix cy (2); 00101 nr_double_t i = 2 * (Id + 2 * Is) * QoverkB / T0 + // shot noise 00102 Kf * qucs::pow (fabs (Id), Af) / qucs::pow (frequency, Ffe) / kB / T0; // flicker noise 00103 cy.set (NODE_C, NODE_C, +i); cy.set (NODE_A, NODE_A, +i); 00104 cy.set (NODE_A, NODE_C, -i); cy.set (NODE_C, NODE_A, -i); 00105 return cy; 00106 } 00107 00108 // Initializes the diode model including temperature and area effects. 00109 void diode::initModel (void) { 00110 // fetch necessary device properties 00111 nr_double_t T = getPropertyDouble ("Temp"); 00112 nr_double_t Tn = getPropertyDouble ("Tnom"); 00113 nr_double_t A = getPropertyDouble ("Area"); 00114 00115 // compute Is temperature and area dependency 00116 nr_double_t Is = getPropertyDouble ("Is"); 00117 nr_double_t N = getPropertyDouble ("N"); 00118 nr_double_t Xti = getPropertyDouble ("Xti"); 00119 nr_double_t Eg = getPropertyDouble ("Eg"); 00120 nr_double_t T1, T2; 00121 T2 = celsius2kelvin (T); 00122 T1 = celsius2kelvin (Tn); 00123 Is = pnCurrent_T (T1, T2, Is, Eg, N, Xti); 00124 setScaledProperty ("Is", Is * A); 00125 00126 // compute Isr temperature and area dependency 00127 nr_double_t Isr = getPropertyDouble ("Isr"); 00128 nr_double_t Nr = getPropertyDouble ("Nr"); 00129 Isr = pnCurrent_T (T1, T2, Isr, Eg, Nr, Xti); 00130 setScaledProperty ("Isr", Isr * A); 00131 00132 // check unphysical parameters 00133 if (Nr < 1.0) { 00134 logprint (LOG_ERROR, "WARNING: Unphysical model parameter Nr = %g in " 00135 "diode `%s'\n", Nr, getName ()); 00136 } 00137 if (N < 1.0) { 00138 logprint (LOG_ERROR, "WARNING: Unphysical model parameter N = %g in " 00139 "diode `%s'\n", N, getName ()); 00140 } 00141 00142 // compute Vj temperature dependency 00143 nr_double_t Vj = getPropertyDouble ("Vj"); 00144 nr_double_t VjT; 00145 VjT = pnPotential_T (T1,T2, Vj); 00146 setScaledProperty ("Vj", VjT); 00147 00148 // compute Cj0 temperature and area dependency 00149 nr_double_t Cj0 = getPropertyDouble ("Cj0"); 00150 nr_double_t M = getPropertyDouble ("M"); 00151 Cj0 = pnCapacitance_T (T1, T2, M, VjT / Vj, Cj0); 00152 setScaledProperty ("Cj0", Cj0 * A); 00153 00154 // check unphysical parameters 00155 if (M > 1.0) { 00156 logprint (LOG_ERROR, "WARNING: Unphysical model parameter M = %g in " 00157 "Diode `%s'\n", M, getName ()); 00158 } 00159 00160 // compute Bv temperature dependency 00161 nr_double_t Bv = getPropertyDouble ("Bv"); 00162 nr_double_t Tbv = getPropertyDouble ("Tbv"); 00163 nr_double_t DT = T2 - T1; 00164 Bv = Bv - Tbv * DT; 00165 setScaledProperty ("Bv", Bv); 00166 00167 // compute Tt temperature dependency 00168 nr_double_t Tt = getPropertyDouble ("Tt"); 00169 nr_double_t Ttt1 = getPropertyDouble ("Ttt1"); 00170 nr_double_t Ttt2 = getPropertyDouble ("Ttt2"); 00171 Tt = Tt * (1 + Ttt1 * DT + Ttt2 * DT * DT); 00172 setScaledProperty ("Tt", Tt); 00173 00174 // compute M temperature dependency 00175 nr_double_t Tm1 = getPropertyDouble ("Tm1"); 00176 nr_double_t Tm2 = getPropertyDouble ("Tm2"); 00177 M = M * (1 + Tm1 * DT + Tm2 * DT * DT); 00178 setScaledProperty ("M", M); 00179 00180 // compute Rs temperature and area dependency 00181 nr_double_t Rs = getPropertyDouble ("Rs"); 00182 nr_double_t Trs = getPropertyDouble ("Trs"); 00183 Rs = Rs * (1 + Trs * DT); 00184 setScaledProperty ("Rs", Rs / A); 00185 } 00186 00187 // Prepares DC (i.e. HB) analysis. 00188 void diode::prepareDC (void) { 00189 // allocate MNA matrices 00190 allocMatrixMNA (); 00191 00192 // initialize scalability 00193 initModel (); 00194 00195 // initialize starting values 00196 Ud = real (getV (NODE_A) - getV (NODE_C)); 00197 for (int i = 0; i < deviceStates (); i++) { 00198 deviceState (i); 00199 UdPrev = Ud; 00200 } 00201 00202 // get device temperature 00203 nr_double_t T = getPropertyDouble ("Temp"); 00204 00205 // possibly insert series resistance 00206 nr_double_t Rs = getScaledProperty ("Rs"); 00207 if (Rs != 0.0) { 00208 // create additional circuit if necessary and reassign nodes 00209 rs = splitResistor (this, rs, "Rs", "anode", NODE_A); 00210 rs->setProperty ("Temp", T); 00211 rs->setProperty ("R", Rs); 00212 rs->setProperty ("Controlled", getName ()); 00213 rs->initDC (); 00214 } 00215 // no series resistance 00216 else { 00217 disableResistor (this, rs, NODE_A); 00218 } 00219 00220 // calculate actual breakdown voltage 00221 Bv = getScaledProperty ("Bv"); 00222 if (Bv != 0) { 00223 nr_double_t Ibv, Is, tol, Ut, Xbv, Xibv; 00224 Ibv = getPropertyDouble ("Ibv"); 00225 Is = getScaledProperty ("Is"); 00226 Ut = celsius2kelvin (T) * kBoverQ; 00227 // adjust very small breakdown currents 00228 if (Ibv < Is * Bv / Ut) { 00229 Ibv = Is * Bv / Ut; 00230 Xbv = Bv; 00231 logprint (LOG_ERROR, "WARNING: Increased breakdown current to %g to " 00232 "match the saturation current %g\n", Ibv, Is); 00233 } 00234 // fit reverse and forward regions 00235 else { 00236 int good = 0; 00237 tol = 1e-3 * Ibv; 00238 Xbv = Bv - Ut * qucs::log (1 + Ibv / Is); 00239 for (int i = 0; i < 25 ; i++) { 00240 Xbv = Bv - Ut * qucs::log (Ibv / Is + 1 - Xbv / Ut); 00241 Xibv = Is * (qucs::exp ((Bv - Xbv) / Ut) - 1 + Xbv / Ut); 00242 if (fabs (Xibv - Ibv) < tol) { 00243 Bv = Xbv; 00244 good = 1; 00245 break; 00246 } 00247 } 00248 if (!good) { 00249 logprint (LOG_ERROR, "WARNING: Unable to fit reverse and forward " 00250 "diode regions using Bv=%g and Ibv=%g\n", Bv, Ibv); 00251 } 00252 } 00253 } 00254 } 00255 00256 // Callback for initializing the DC analysis. 00257 void diode::initDC (void) { 00258 deviceStates (StateVars, 1); 00259 doHB = false; 00260 prepareDC (); 00261 } 00262 00263 // Callback for restarting the DC analysis. 00264 void diode::restartDC (void) { 00265 // apply starting value to previous iteration value 00266 UdPrev = real (getV (NODE_A) - getV (NODE_C)); 00267 } 00268 00269 // Callback for DC analysis. 00270 void diode::calcDC (void) { 00271 // get device properties 00272 nr_double_t Is = getScaledProperty ("Is"); 00273 nr_double_t N = getPropertyDouble ("N"); 00274 nr_double_t Isr = getScaledProperty ("Isr"); 00275 nr_double_t Nr = getPropertyDouble ("Nr"); 00276 nr_double_t Ikf = getPropertyDouble ("Ikf"); 00277 nr_double_t T = getPropertyDouble ("Temp"); 00278 00279 nr_double_t Ut, Ieq, Ucrit, gtiny; 00280 00281 T = celsius2kelvin (T); 00282 Ut = T * kBoverQ; 00283 Ud = real (getV (NODE_A) - getV (NODE_C)); 00284 00285 // critical voltage necessary for bad start values 00286 Ucrit = pnCriticalVoltage (Is, N * Ut); 00287 if (Bv != 0 && Ud < std::min (0.0, -Bv + 10 * N * Ut)) { 00288 nr_double_t V = -(Ud + Bv); 00289 V = pnVoltage (V, -(UdPrev + Bv), Ut * N, Ucrit); 00290 Ud = -(V + Bv); 00291 } 00292 else { 00293 Ud = pnVoltage (Ud, UdPrev, Ut * N, Ucrit); 00294 } 00295 UdPrev = Ud; 00296 00297 // tiny derivative for little junction voltage 00298 gtiny = (Ud < - 10 * Ut * N && Bv != 0) ? (Is + Isr) : 0; 00299 00300 if (Ud >= -3 * N * Ut) { // forward region 00301 gd = pnConductance (Ud, Is, Ut * N) + pnConductance (Ud, Isr, Ut * Nr); 00302 Id = pnCurrent (Ud, Is, Ut * N) + pnCurrent (Ud, Isr, Ut * Nr); 00303 } 00304 else if (Bv == 0 || Ud >= -Bv) { // reverse region 00305 nr_double_t a = 3 * N * Ut / (Ud * euler); 00306 a = cubic (a); 00307 Id = -Is * (1 + a); 00308 gd = +Is * 3 * a / Ud; 00309 } 00310 else { // middle region 00311 nr_double_t a = qucs::exp (-(Bv + Ud) / N / Ut); 00312 Id = -Is * a; 00313 gd = +Is * a / Ut / N; 00314 } 00315 00316 // knee current calculations 00317 if (Ikf != 0.0) { 00318 nr_double_t a = Ikf / (Ikf + Id); 00319 gd *= 0.5 * (2 - Id * a / Ikf) * qucs::sqrt (a); 00320 Id *= qucs::sqrt (a); 00321 } 00322 00323 Id += gtiny * Ud; 00324 gd += gtiny; 00325 00326 // HB simulation 00327 if (doHB) { 00328 Ieq = Id; 00329 setGV (NODE_C, -gd * Ud); 00330 setGV (NODE_A, +gd * Ud); 00331 } 00332 // DC and transient simulation 00333 else { 00334 Ieq = Id - Ud * gd; 00335 } 00336 00337 // fill in I-Vector 00338 setI (NODE_C, +Ieq); 00339 setI (NODE_A, -Ieq); 00340 00341 // fill in G-Matrix 00342 setY (NODE_C, NODE_C, +gd); setY (NODE_A, NODE_A, +gd); 00343 setY (NODE_C, NODE_A, -gd); setY (NODE_A, NODE_C, -gd); 00344 } 00345 00346 // Saves operating points (voltages). 00347 void diode::saveOperatingPoints (void) { 00348 nr_double_t Vd = real (getV (NODE_A) - getV (NODE_C)); 00349 setOperatingPoint ("Vd", Vd); 00350 } 00351 00352 // Loads operating points (voltages). 00353 void diode::loadOperatingPoints (void) { 00354 Ud = getOperatingPoint ("Vd"); 00355 } 00356 00357 // Calculates and saves operating points. 00358 void diode::calcOperatingPoints (void) { 00359 00360 // load operating points 00361 loadOperatingPoints (); 00362 00363 // get necessary properties 00364 nr_double_t M = getScaledProperty ("M"); 00365 nr_double_t Cj0 = getScaledProperty ("Cj0"); 00366 nr_double_t Vj = getScaledProperty ("Vj"); 00367 nr_double_t Fc = getPropertyDouble ("Fc"); 00368 nr_double_t Cp = getPropertyDouble ("Cp"); 00369 nr_double_t Tt = getScaledProperty ("Tt"); 00370 00371 // calculate capacitances and charges 00372 nr_double_t Cd; 00373 Cd = pnCapacitance (Ud, Cj0, Vj, M, Fc) + Tt * gd + Cp; 00374 Qd = pnCharge (Ud, Cj0, Vj, M, Fc) + Tt * Id + Cp * Ud; 00375 00376 // save operating points 00377 setOperatingPoint ("gd", gd); 00378 setOperatingPoint ("Id", Id); 00379 setOperatingPoint ("Cd", Cd); 00380 } 00381 00382 // Callback for initializing the AC analysis. 00383 void diode::initAC (void) { 00384 allocMatrixMNA (); 00385 } 00386 00387 // Callback for the AC analysis. 00388 void diode::calcAC (nr_double_t frequency) { 00389 nr_double_t gd = getOperatingPoint ("gd"); 00390 nr_double_t Cd = getOperatingPoint ("Cd"); 00391 nr_complex_t y = nr_complex_t (gd, Cd * 2.0 * pi * frequency); 00392 setY (NODE_C, NODE_C, +y); setY (NODE_A, NODE_A, +y); 00393 setY (NODE_C, NODE_A, -y); setY (NODE_A, NODE_C, -y); 00394 } 00395 00396 // Callback for the AC noise analysis. 00397 void diode::calcNoiseAC (nr_double_t frequency) { 00398 setMatrixN (calcMatrixCy (frequency)); 00399 } 00400 00401 #define qState 0 // charge state 00402 #define cState 1 // current state 00403 00404 // Callback for initializing the TR analysis. 00405 void diode::initTR (void) { 00406 setStates (2); 00407 initDC (); 00408 } 00409 00410 // Callback for the TR analysis. 00411 void diode::calcTR (nr_double_t) { 00412 calcDC (); 00413 saveOperatingPoints (); 00414 calcOperatingPoints (); 00415 00416 nr_double_t Cd = getOperatingPoint ("Cd"); 00417 00418 transientCapacitance (qState, NODE_A, NODE_C, Cd, Ud, Qd); 00419 } 00420 00421 // Callback for initializing the HB analysis. 00422 void diode::initHB (int frequencies) { 00423 deviceStates (StateVars, frequencies); 00424 doHB = true; 00425 prepareDC (); 00426 allocMatrixHB (); 00427 } 00428 00429 // Callback for the HB analysis. 00430 void diode::calcHB (int frequency) { 00431 // set current frequency state 00432 deviceState (frequency); 00433 00434 // g's (dI/dU) into Y-Matrix and I's into I-Vector 00435 calcDC (); 00436 00437 // calculate Q and C 00438 saveOperatingPoints (); 00439 calcOperatingPoints (); 00440 00441 nr_double_t Cd = getOperatingPoint ("Cd"); 00442 00443 // fill in Q's in Q-Vector 00444 setQ (NODE_C, +Qd); 00445 setQ (NODE_A, -Qd); 00446 00447 setCV (NODE_C, -Cd * Ud); 00448 setCV (NODE_A, +Cd * Ud); 00449 00450 // fill in C's (dQ/dU) into QV-Matrix 00451 setQV (NODE_C, NODE_C, +Cd); setQV (NODE_A, NODE_A, +Cd); 00452 setQV (NODE_C, NODE_A, -Cd); setQV (NODE_A, NODE_C, -Cd); 00453 } 00454 00455 // properties 00456 PROP_REQ [] = { 00457 { "Is", PROP_REAL, { 1e-15, PROP_NO_STR }, PROP_POS_RANGE }, 00458 { "N", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1e-6, 100) }, 00459 { "M", PROP_REAL, { 0.5, PROP_NO_STR }, PROP_RNGII (0, 2) }, 00460 { "Cj0", PROP_REAL, { 10e-15, PROP_NO_STR }, PROP_POS_RANGE }, 00461 { "Vj", PROP_REAL, { 0.7, PROP_NO_STR }, PROP_RNGXI (0, 10) }, 00462 PROP_NO_PROP }; 00463 PROP_OPT [] = { 00464 { "Rs", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00465 { "Isr", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00466 { "Nr", PROP_REAL, { 2, PROP_NO_STR }, PROP_RNGII (0.1, 100) }, 00467 { "Bv", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00468 { "Ibv", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE }, 00469 { "Ikf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00470 { "Tt", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00471 { "Fc", PROP_REAL, { 0.5, PROP_NO_STR }, PROP_RNGIX (0, 1) }, 00472 { "Cp", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00473 { "Kf", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00474 { "Af", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE }, 00475 { "Ffe", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGE }, 00476 { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) }, 00477 { "Xti", PROP_REAL, { 3, PROP_NO_STR }, PROP_POS_RANGE }, 00478 { "Eg", PROP_REAL, { EgSi, PROP_NO_STR }, PROP_POS_RANGE }, 00479 { "Tbv", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE }, 00480 { "Trs", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00481 { "Ttt1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00482 { "Ttt2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00483 { "Tm1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00484 { "Tm2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00485 { "Tnom", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) }, 00486 { "Area", PROP_REAL, { 1, PROP_NO_STR }, PROP_POS_RANGEX }, 00487 PROP_NO_PROP }; 00488 struct define_t diode::cirdef = 00489 { "Diode", 2, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_NONLINEAR, PROP_DEF };