1 /********** 2 Copyright 1991 Regents of the University of California. All rights reserved. 3 Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group 4 Author: 1991 David A. Gates, U. C. Berkeley CAD Group 5 **********/ 6 7 #include "ngspice/ngspice.h" 8 #include "ngspice/numglobs.h" 9 #include "ngspice/numenum.h" 10 #include "ngspice/twomesh.h" 11 #include "ngspice/twodev.h" 12 #include "ngspice/bool.h" 13 #include "ngspice/spmatrix.h" 14 #include "twoddefs.h" 15 #include "twodext.h" 16 #include "ngspice/cidersupt.h" 17 #include "../../maths/misc/bernoull.h" 18 19 /* 20 * Functions to setup and solve the continuity equations. 21 * Both continuity equations are solved. 22 * Separate functions are used for one continuity equation. 23 */ 24 25 26 /* 27 * setup matrix pointers to Jacobian values and 28 * store direct pointers with the nodes 29 */ 30 31 void TWO_jacBuild(TWOdevice * pDevice)32 TWO_jacBuild(TWOdevice *pDevice) 33 { 34 SMPmatrix *matrix = pDevice->matrix; 35 TWOelem *pElem; 36 TWOnode *pNode; 37 TWOchannel *pCh; 38 int eIndex, nIndex; 39 int nextIndex; /* index of node to find next element */ 40 int psiEqn, nEqn, pEqn; /* scratch for deref'd eqn numbers */ 41 int psiEqnTL = 0, nEqnTL = 0, pEqnTL = 0; 42 int psiEqnTR = 0, nEqnTR = 0, pEqnTR = 0; 43 int psiEqnBR = 0, nEqnBR = 0, pEqnBR = 0; 44 int psiEqnBL = 0, nEqnBL = 0, pEqnBL = 0; 45 int psiEqnInM = 0, psiEqnInP = 0; /* scratch for deref'd surface eqns */ 46 int psiEqnOxM = 0, psiEqnOxP = 0; /* M= more negative, P= more positive */ 47 48 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 49 pElem = pDevice->elements[ eIndex ]; 50 51 /* first the self terms */ 52 for ( nIndex = 0; nIndex <= 3; nIndex++ ) { 53 pNode = pElem->pNodes[ nIndex ]; 54 /* get poisson-only pointer */ 55 psiEqn = pNode->psiEqn; 56 pNode->fPsiPsi = spGetElement( matrix, psiEqn, psiEqn ); 57 58 if ( pElem->elemType == SEMICON ) { 59 /* get continuity-coupling terms */ 60 nEqn = pNode->nEqn; 61 pEqn = pNode->pEqn; 62 /* pointers for additional terms */ 63 pNode->fPsiN = spGetElement( matrix, psiEqn, nEqn ); 64 pNode->fPsiP = spGetElement( matrix, psiEqn, pEqn ); 65 pNode->fNPsi = spGetElement( matrix, nEqn, psiEqn ); 66 pNode->fNN = spGetElement( matrix, nEqn, nEqn ); 67 pNode->fNP = spGetElement( matrix, nEqn, pEqn ); 68 pNode->fPPsi = spGetElement( matrix, pEqn, psiEqn ); 69 pNode->fPN = spGetElement( matrix, pEqn, nEqn ); 70 pNode->fPP = spGetElement( matrix, pEqn, pEqn ); 71 } else { 72 nEqn = 0; 73 pEqn = 0; 74 } 75 /* save equation indices */ 76 switch ( nIndex ) { 77 case 0: /* TL Node */ 78 psiEqnTL = psiEqn; 79 nEqnTL = nEqn; 80 pEqnTL = pEqn; 81 break; 82 case 1: /* TR Node */ 83 psiEqnTR = psiEqn; 84 nEqnTR = nEqn; 85 pEqnTR = pEqn; 86 break; 87 case 2: /* BR Node */ 88 psiEqnBR = psiEqn; 89 nEqnBR = nEqn; 90 pEqnBR = pEqn; 91 break; 92 case 3: /* BL Node */ 93 psiEqnBL = psiEqn; 94 nEqnBL = nEqn; 95 pEqnBL = pEqn; 96 break; 97 default: 98 break; 99 } 100 } 101 102 /* now terms to couple to adjacent nodes */ 103 pNode = pElem->pTLNode; 104 pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnTL, psiEqnTR ); 105 pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTL, psiEqnBL ); 106 if ( pElem->elemType == SEMICON ) { 107 /* continuity equation pointers */ 108 pNode->fNPsiiP1 = spGetElement( matrix, nEqnTL, psiEqnTR ); 109 pNode->fNNiP1 = spGetElement( matrix, nEqnTL, nEqnTR ); 110 pNode->fNPsijP1 = spGetElement( matrix, nEqnTL, psiEqnBL ); 111 pNode->fNNjP1 = spGetElement( matrix, nEqnTL, nEqnBL ); 112 pNode->fPPsiiP1 = spGetElement( matrix, pEqnTL, psiEqnTR ); 113 pNode->fPPiP1 = spGetElement( matrix, pEqnTL, pEqnTR ); 114 pNode->fPPsijP1 = spGetElement( matrix, pEqnTL, psiEqnBL ); 115 pNode->fPPjP1 = spGetElement( matrix, pEqnTL, pEqnBL ); 116 /* Surface Mobility Model depends on diagonal node values */ 117 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 118 pNode->fNPsiiP1jP1 = spGetElement( matrix, nEqnTL, psiEqnBR ); 119 pNode->fNNiP1jP1 = spGetElement( matrix, nEqnTL, nEqnBR ); 120 pNode->fPPsiiP1jP1 = spGetElement( matrix, pEqnTL, psiEqnBR ); 121 pNode->fPPiP1jP1 = spGetElement( matrix, pEqnTL, pEqnBR ); 122 } 123 } 124 125 pNode = pElem->pTRNode; 126 pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnTR, psiEqnTL ); 127 pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTR, psiEqnBR ); 128 if ( pElem->elemType == SEMICON ) { 129 /* continuity equation pointers */ 130 pNode->fNPsiiM1 = spGetElement( matrix, nEqnTR, psiEqnTL ); 131 pNode->fNNiM1 = spGetElement( matrix, nEqnTR, nEqnTL ); 132 pNode->fNPsijP1 = spGetElement( matrix, nEqnTR, psiEqnBR ); 133 pNode->fNNjP1 = spGetElement( matrix, nEqnTR, nEqnBR ); 134 pNode->fPPsiiM1 = spGetElement( matrix, pEqnTR, psiEqnTL ); 135 pNode->fPPiM1 = spGetElement( matrix, pEqnTR, pEqnTL ); 136 pNode->fPPsijP1 = spGetElement( matrix, pEqnTR, psiEqnBR ); 137 pNode->fPPjP1 = spGetElement( matrix, pEqnTR, pEqnBR ); 138 /* Surface Mobility Model depends on diagonal node values */ 139 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 140 pNode->fNPsiiM1jP1 = spGetElement( matrix, nEqnTR, psiEqnBL ); 141 pNode->fNNiM1jP1 = spGetElement( matrix, nEqnTR, nEqnBL ); 142 pNode->fPPsiiM1jP1 = spGetElement( matrix, pEqnTR, psiEqnBL ); 143 pNode->fPPiM1jP1 = spGetElement( matrix, pEqnTR, pEqnBL ); 144 } 145 } 146 147 pNode = pElem->pBRNode; 148 pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnBR, psiEqnBL ); 149 pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBR, psiEqnTR ); 150 if ( pElem->elemType == SEMICON ) { 151 /* continuity equation pointers */ 152 pNode->fNPsiiM1 = spGetElement( matrix, nEqnBR, psiEqnBL ); 153 pNode->fNNiM1 = spGetElement( matrix, nEqnBR, nEqnBL ); 154 pNode->fNPsijM1 = spGetElement( matrix, nEqnBR, psiEqnTR ); 155 pNode->fNNjM1 = spGetElement( matrix, nEqnBR, nEqnTR ); 156 pNode->fPPsiiM1 = spGetElement( matrix, pEqnBR, psiEqnBL ); 157 pNode->fPPiM1 = spGetElement( matrix, pEqnBR, pEqnBL ); 158 pNode->fPPsijM1 = spGetElement( matrix, pEqnBR, psiEqnTR ); 159 pNode->fPPjM1 = spGetElement( matrix, pEqnBR, pEqnTR ); 160 /* Surface Mobility Model depends on diagonal node values */ 161 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 162 pNode->fNPsiiM1jM1 = spGetElement( matrix, nEqnBR, psiEqnTL ); 163 pNode->fNNiM1jM1 = spGetElement( matrix, nEqnBR, nEqnTL ); 164 pNode->fPPsiiM1jM1 = spGetElement( matrix, pEqnBR, psiEqnTL ); 165 pNode->fPPiM1jM1 = spGetElement( matrix, pEqnBR, pEqnTL ); 166 } 167 } 168 169 pNode = pElem->pBLNode; 170 pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnBL, psiEqnBR ); 171 pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBL, psiEqnTL ); 172 if ( pElem->elemType == SEMICON ) { 173 /* continuity equation pointers */ 174 pNode->fNPsiiP1 = spGetElement( matrix, nEqnBL, psiEqnBR ); 175 pNode->fNNiP1 = spGetElement( matrix, nEqnBL, nEqnBR ); 176 pNode->fNPsijM1 = spGetElement( matrix, nEqnBL, psiEqnTL ); 177 pNode->fNNjM1 = spGetElement( matrix, nEqnBL, nEqnTL ); 178 pNode->fPPsiiP1 = spGetElement( matrix, pEqnBL, psiEqnBR ); 179 pNode->fPPiP1 = spGetElement( matrix, pEqnBL, pEqnBR ); 180 pNode->fPPsijM1 = spGetElement( matrix, pEqnBL, psiEqnTL ); 181 pNode->fPPjM1 = spGetElement( matrix, pEqnBL, pEqnTL ); 182 /* Surface Mobility Model depends on diagonal node values */ 183 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 184 pNode->fNPsiiP1jM1 = spGetElement( matrix, nEqnBL, psiEqnTR ); 185 pNode->fNNiP1jM1 = spGetElement( matrix, nEqnBL, nEqnTR ); 186 pNode->fPPsiiP1jM1 = spGetElement( matrix, pEqnBL, psiEqnTR ); 187 pNode->fPPiP1jM1 = spGetElement( matrix, pEqnBL, pEqnTR ); 188 } 189 } 190 } 191 /* 192 * Add terms for surface-field of inversion-layer mobility model. 193 * Elements MUST be made from silicon for this to work. 194 * No empty elements are allowed. 195 * Don't need these pointers if SurfaceMobility isn't set. 196 */ 197 if ( MobDeriv && SurfaceMobility ) { 198 for ( pCh = pDevice->pChannel; pCh != NULL; 199 pCh = pCh->next ) { 200 pElem = pCh->pNElem; 201 switch (pCh->type) { 202 case 0: 203 psiEqnInM = pElem->pBLNode->psiEqn; 204 psiEqnInP = pElem->pBRNode->psiEqn; 205 psiEqnOxM = pElem->pTLNode->psiEqn; 206 psiEqnOxP = pElem->pTRNode->psiEqn; 207 break; 208 case 1: 209 psiEqnInM = pElem->pTLNode->psiEqn; 210 psiEqnInP = pElem->pBLNode->psiEqn; 211 psiEqnOxM = pElem->pTRNode->psiEqn; 212 psiEqnOxP = pElem->pBRNode->psiEqn; 213 break; 214 case 2: 215 psiEqnInM = pElem->pTLNode->psiEqn; 216 psiEqnInP = pElem->pTRNode->psiEqn; 217 psiEqnOxM = pElem->pBLNode->psiEqn; 218 psiEqnOxP = pElem->pBRNode->psiEqn; 219 break; 220 case 3: 221 psiEqnInM = pElem->pTRNode->psiEqn; 222 psiEqnInP = pElem->pBRNode->psiEqn; 223 psiEqnOxM = pElem->pTLNode->psiEqn; 224 psiEqnOxP = pElem->pBLNode->psiEqn; 225 break; 226 } 227 pElem = pCh->pSeed; 228 nextIndex = (pCh->type + 2)%4; 229 while (pElem && pElem->channel == pCh->id) { 230 for ( nIndex = 0; nIndex <= 3; nIndex++ ) { 231 pNode = pElem->pNodes[ nIndex ]; 232 psiEqn = pNode->psiEqn; 233 nEqn = pNode->nEqn; 234 pEqn = pNode->pEqn; 235 if ( pCh->type % 2 == 0 ) { /* Vertical Slice */ 236 if ( nIndex == 0 || nIndex == 3 ) { /* Left Side */ 237 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInM ); 238 pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP ); 239 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxM ); 240 pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP ); 241 pNode->fPPsiIn = spGetElement( matrix, pEqn, psiEqnInM ); 242 pNode->fPPsiInP1 = spGetElement( matrix, pEqn, psiEqnInP ); 243 pNode->fPPsiOx = spGetElement( matrix, pEqn, psiEqnOxM ); 244 pNode->fPPsiOxP1 = spGetElement( matrix, pEqn, psiEqnOxP ); 245 } else { /* Right Side */ 246 pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM ); 247 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInP ); 248 pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM ); 249 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxP ); 250 pNode->fPPsiInM1 = spGetElement( matrix, pEqn, psiEqnInM ); 251 pNode->fPPsiIn = spGetElement( matrix, pEqn, psiEqnInP ); 252 pNode->fPPsiOxM1 = spGetElement( matrix, pEqn, psiEqnOxM ); 253 pNode->fPPsiOx = spGetElement( matrix, pEqn, psiEqnOxP ); 254 } 255 } else { /* Horizontal Slice */ 256 if ( nIndex == 0 || nIndex == 3 ) { /* Left (Top?) Side : bug 483 */ 257 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInM ); 258 pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP ); 259 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxM ); 260 pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP ); 261 pNode->fPPsiIn = spGetElement( matrix, pEqn, psiEqnInM ); 262 pNode->fPPsiInP1 = spGetElement( matrix, pEqn, psiEqnInP ); 263 pNode->fPPsiOx = spGetElement( matrix, pEqn, psiEqnOxM ); 264 pNode->fPPsiOxP1 = spGetElement( matrix, pEqn, psiEqnOxP ); 265 } else { /* Bottom Side */ 266 pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM ); 267 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInP ); 268 pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM ); 269 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxP ); 270 pNode->fPPsiInM1 = spGetElement( matrix, pEqn, psiEqnInM ); 271 pNode->fPPsiIn = spGetElement( matrix, pEqn, psiEqnInP ); 272 pNode->fPPsiOxM1 = spGetElement( matrix, pEqn, psiEqnOxM ); 273 pNode->fPPsiOx = spGetElement( matrix, pEqn, psiEqnOxP ); 274 } 275 } 276 } /* endfor nIndex */ 277 pElem = pElem->pElems[ nextIndex ]; 278 } /* endwhile pElem */ 279 } /* endfor pCh */ 280 } /* endif SurfaceMobility */ 281 } 282 283 284 /* 285 * The Jacobian and Rhs are loaded by the following function. 286 * Inputs are the transient analysis flag and the transient 287 * information structure 288 */ 289 290 void TWO_sysLoad(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)291 TWO_sysLoad(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info) 292 { 293 TWOelem *pElem; 294 TWOnode *pNode; 295 TWOedge *pHEdge, *pVEdge; 296 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 297 TWOchannel *pCh; 298 int index, eIndex; 299 int nextIndex; /* index of node to find next element */ 300 double *pRhs = pDevice->rhs; 301 double dx, dy, dxdy, dyOverDx, dxOverDy; 302 double ds; 303 double dPsiT, dPsiB, dPsiL, dPsiR; 304 double rhsN, rhsP; 305 double generation; 306 double nConc, pConc; 307 double perTime = 0.0; 308 309 /* first compute the currents and derivatives */ 310 TWO_commonTerms( pDevice, FALSE, tranAnalysis, info ); 311 312 /* find reciprocal timestep */ 313 if ( tranAnalysis ) { 314 perTime = info->intCoeff[0]; 315 } 316 317 /* zero the rhs vector */ 318 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) { 319 pRhs[ index ] = 0.0; 320 } 321 322 /* zero the matrix */ 323 spClear( pDevice->matrix ); 324 325 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 326 pElem = pDevice->elements[ eIndex ]; 327 328 dx = 0.5 * pElem->dx; 329 dy = 0.5 * pElem->dy; 330 dxdy = dx * dy; 331 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 332 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 333 334 pTEdge = pElem->pTopEdge; 335 pBEdge = pElem->pBotEdge; 336 pLEdge = pElem->pLeftEdge; 337 pREdge = pElem->pRightEdge; 338 dPsiT = pTEdge->dPsi; 339 dPsiB = pBEdge->dPsi; 340 dPsiL = pLEdge->dPsi; 341 dPsiR = pREdge->dPsi; 342 343 /* load for all i,j */ 344 for ( index = 0; index <= 3; index++ ) { 345 pNode = pElem->pNodes[ index ]; 346 if ( pNode->nodeType != CONTACT ) { 347 *(pNode->fPsiPsi) += dyOverDx + dxOverDy; 348 if ( index <= 1 ) { 349 pHEdge = pTEdge; 350 } else { 351 pHEdge = pBEdge; 352 } 353 if ( index == 0 || index == 3 ) { 354 pVEdge = pLEdge; 355 } else { 356 pVEdge = pREdge; 357 } 358 /* Add surface state charges. */ 359 pRhs[ pNode->psiEqn ] += dx * pHEdge->qf; 360 pRhs[ pNode->psiEqn ] += dy * pVEdge->qf; 361 if ( pElem->elemType == SEMICON ) { 362 *(pNode->fPsiN) += dxdy; 363 *(pNode->fPsiP) -= dxdy; 364 *(pNode->fNPsi) -= dy * pHEdge->dJnDpsiP1 + 365 dx * pVEdge->dJnDpsiP1; 366 *(pNode->fPPsi) -= dy * pHEdge->dJpDpsiP1 + 367 dx * pVEdge->dJpDpsiP1; 368 369 nConc = pDevice->devState0 [pNode->nodeN]; 370 pConc = pDevice->devState0 [pNode->nodeP]; 371 pRhs[ pNode->psiEqn ] += dxdy * (pNode->netConc + pConc - nConc); 372 373 /* Handle generation terms */ 374 *(pNode->fNN) -= dxdy * pNode->dUdN; 375 *(pNode->fNP) -= dxdy * pNode->dUdP; 376 *(pNode->fPP) += dxdy * pNode->dUdP; 377 *(pNode->fPN) += dxdy * pNode->dUdN; 378 rhsN = - dxdy * pNode->uNet; 379 rhsP = dxdy * pNode->uNet; 380 if ( AvalancheGen ) { 381 generation = TWOavalanche( pElem, pNode ); 382 rhsN += dxdy * generation; 383 rhsP -= dxdy * generation; 384 } 385 pRhs[ pNode->nEqn ] -= rhsN; 386 pRhs[ pNode->pEqn ] -= rhsP; 387 388 /* Handle dXdT continuity terms */ 389 if ( tranAnalysis ) { 390 *(pNode->fNN) -= dxdy * perTime; 391 *(pNode->fPP) += dxdy * perTime; 392 pRhs[ pNode->nEqn ] += dxdy * pNode->dNdT; 393 pRhs[ pNode->pEqn ] -= dxdy * pNode->dPdT; 394 } 395 } 396 } 397 } 398 399 /* Handle neighbor and edge dependent terms */ 400 pNode = pElem->pTLNode; 401 if ( pNode->nodeType != CONTACT ) { 402 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL; 403 *(pNode->fPsiPsiiP1) -= dyOverDx; 404 *(pNode->fPsiPsijP1) -= dxOverDy; 405 if ( pElem->elemType == SEMICON ) { 406 pRhs[ pNode->nEqn ] -= dy * pTEdge->jn + dx * pLEdge->jn; 407 pRhs[ pNode->pEqn ] -= dy * pTEdge->jp + dx * pLEdge->jp; 408 *(pNode->fNN) += dy * pTEdge->dJnDn + dx * pLEdge->dJnDn; 409 *(pNode->fPP) += dy * pTEdge->dJpDp + dx * pLEdge->dJpDp; 410 *(pNode->fNPsiiP1) += dy * pTEdge->dJnDpsiP1; 411 *(pNode->fNNiP1) += dy * pTEdge->dJnDnP1; 412 *(pNode->fPPsiiP1) += dy * pTEdge->dJpDpsiP1; 413 *(pNode->fPPiP1) += dy * pTEdge->dJpDpP1; 414 *(pNode->fNPsijP1) += dx * pLEdge->dJnDpsiP1; 415 *(pNode->fNNjP1) += dx * pLEdge->dJnDnP1; 416 *(pNode->fPPsijP1) += dx * pLEdge->dJpDpsiP1; 417 *(pNode->fPPjP1) += dx * pLEdge->dJpDpP1; 418 } 419 } 420 pNode = pElem->pTRNode; 421 if ( pNode->nodeType != CONTACT ) { 422 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR; 423 *(pNode->fPsiPsiiM1) -= dyOverDx; 424 *(pNode->fPsiPsijP1) -= dxOverDy; 425 if ( pElem->elemType == SEMICON ) { 426 pRhs[ pNode->nEqn ] -= -dy * pTEdge->jn + dx * pREdge->jn; 427 pRhs[ pNode->pEqn ] -= -dy * pTEdge->jp + dx * pREdge->jp; 428 *(pNode->fNN) += -dy * pTEdge->dJnDnP1 + dx * pREdge->dJnDn; 429 *(pNode->fPP) += -dy * pTEdge->dJpDpP1 + dx * pREdge->dJpDp; 430 *(pNode->fNPsiiM1) += dy * pTEdge->dJnDpsiP1; 431 *(pNode->fNNiM1) -= dy * pTEdge->dJnDn; 432 *(pNode->fPPsiiM1) += dy * pTEdge->dJpDpsiP1; 433 *(pNode->fPPiM1) -= dy * pTEdge->dJpDp; 434 *(pNode->fNPsijP1) += dx * pREdge->dJnDpsiP1; 435 *(pNode->fNNjP1) += dx * pREdge->dJnDnP1; 436 *(pNode->fPPsijP1) += dx * pREdge->dJpDpsiP1; 437 *(pNode->fPPjP1) += dx * pREdge->dJpDpP1; 438 } 439 } 440 pNode = pElem->pBRNode; 441 if ( pNode->nodeType != CONTACT ) { 442 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR; 443 *(pNode->fPsiPsiiM1) -= dyOverDx; 444 *(pNode->fPsiPsijM1) -= dxOverDy; 445 if ( pElem->elemType == SEMICON ) { 446 pRhs[ pNode->nEqn ] -= -dy * pBEdge->jn - dx * pREdge->jn; 447 pRhs[ pNode->pEqn ] -= -dy * pBEdge->jp - dx * pREdge->jp; 448 *(pNode->fNN) += -dy * pBEdge->dJnDnP1 - dx * pREdge->dJnDnP1; 449 *(pNode->fPP) += -dy * pBEdge->dJpDpP1 - dx * pREdge->dJpDpP1; 450 *(pNode->fNPsiiM1) += dy * pBEdge->dJnDpsiP1; 451 *(pNode->fNNiM1) -= dy * pBEdge->dJnDn; 452 *(pNode->fPPsiiM1) += dy * pBEdge->dJpDpsiP1; 453 *(pNode->fPPiM1) -= dy * pBEdge->dJpDp; 454 *(pNode->fNPsijM1) += dx * pREdge->dJnDpsiP1; 455 *(pNode->fNNjM1) -= dx * pREdge->dJnDn; 456 *(pNode->fPPsijM1) += dx * pREdge->dJpDpsiP1; 457 *(pNode->fPPjM1) -= dx * pREdge->dJpDp; 458 } 459 } 460 pNode = pElem->pBLNode; 461 if ( pNode->nodeType != CONTACT ) { 462 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL; 463 *(pNode->fPsiPsiiP1) -= dyOverDx; 464 *(pNode->fPsiPsijM1) -= dxOverDy; 465 if ( pElem->elemType == SEMICON ) { 466 pRhs[ pNode->nEqn ] -= dy * pBEdge->jn - dx * pLEdge->jn; 467 pRhs[ pNode->pEqn ] -= dy * pBEdge->jp - dx * pLEdge->jp; 468 *(pNode->fNN) += dy * pBEdge->dJnDn - dx * pLEdge->dJnDnP1; 469 *(pNode->fPP) += dy * pBEdge->dJpDp - dx * pLEdge->dJpDpP1; 470 *(pNode->fNPsiiP1) += dy * pBEdge->dJnDpsiP1; 471 *(pNode->fNNiP1) += dy * pBEdge->dJnDnP1; 472 *(pNode->fPPsiiP1) += dy * pBEdge->dJpDpsiP1; 473 *(pNode->fPPiP1) += dy * pBEdge->dJpDpP1; 474 *(pNode->fNPsijM1) += dx * pLEdge->dJnDpsiP1; 475 *(pNode->fNNjM1) -= dx * pLEdge->dJnDn; 476 *(pNode->fPPsijM1) += dx * pLEdge->dJpDpsiP1; 477 *(pNode->fPPjM1) -= dx * pLEdge->dJpDp; 478 } 479 } 480 } 481 482 /* Calculate the Inversion-Layer Mobility Dependent Terms in Jac. */ 483 if ( MobDeriv && SurfaceMobility ) { 484 for ( pCh = pDevice->pChannel; pCh != NULL; 485 pCh = pCh->next ) { 486 /* Find effective height of oxide element at interface. */ 487 if ( pCh->type%2 == 0 ) { /* Vertical slice */ 488 ds = pCh->pNElem->dy / pCh->pNElem->epsRel; 489 } else { /* Horizontal slice */ 490 ds = pCh->pNElem->dx / pCh->pNElem->epsRel; 491 } 492 pElem = pCh->pSeed; 493 nextIndex = (pCh->type + 2)%4; 494 while (pElem && pElem->channel == pCh->id) { 495 TWO_mobDeriv( pElem, pCh->type, ds ); 496 pElem = pElem->pElems[ nextIndex ]; 497 } 498 } /* endfor pCh != NULL */ 499 } /* endif MobDeriv and SurfaceMobility */ 500 } 501 502 503 /* this function used only for direct method ac analysis 504 Used to load only the dc Jacobian matrix. Rhs is unaffected 505 */ 506 507 void TWO_jacLoad(TWOdevice * pDevice)508 TWO_jacLoad(TWOdevice *pDevice) 509 { 510 TWOelem *pElem; 511 TWOnode *pNode; 512 TWOedge *pHEdge, *pVEdge; 513 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 514 TWOchannel *pCh; 515 int index, eIndex; 516 int nextIndex; /* index of node to find next element */ 517 double dx, dy, dxdy, dyOverDx, dxOverDy; 518 double ds; 519 520 /* first compute the currents and derivatives */ 521 TWO_commonTerms( pDevice, FALSE, FALSE, NULL ); 522 523 /* zero the matrix */ 524 spClear( pDevice->matrix ); 525 526 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 527 pElem = pDevice->elements[ eIndex ]; 528 dx = 0.5 * pElem->dx; 529 dy = 0.5 * pElem->dy; 530 dxdy = dx * dy; 531 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 532 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 533 534 pTEdge = pElem->pTopEdge; 535 pBEdge = pElem->pBotEdge; 536 pLEdge = pElem->pLeftEdge; 537 pREdge = pElem->pRightEdge; 538 539 /* load for all i,j */ 540 for ( index = 0; index <= 3; index++ ) { 541 pNode = pElem->pNodes[ index ]; 542 if ( pNode->nodeType != CONTACT ) { 543 *(pNode->fPsiPsi) += dyOverDx + dxOverDy; 544 if ( pElem->elemType == SEMICON ) { 545 if ( index <= 1 ) { 546 pHEdge = pTEdge; 547 } else { 548 pHEdge = pBEdge; 549 } 550 if ( index == 0 || index == 3 ) { 551 pVEdge = pLEdge; 552 } else { 553 pVEdge = pREdge; 554 } 555 *(pNode->fPsiN) += dxdy; 556 *(pNode->fPsiP) -= dxdy; 557 *(pNode->fNPsi) -= dy * pHEdge->dJnDpsiP1 + 558 dx * pVEdge->dJnDpsiP1; 559 *(pNode->fPPsi) -= dy * pHEdge->dJpDpsiP1 + 560 dx * pVEdge->dJpDpsiP1; 561 562 /* Handle generation terms */ 563 *(pNode->fNN) -= dxdy * pNode->dUdN; 564 *(pNode->fNP) -= dxdy * pNode->dUdP; 565 *(pNode->fPP) += dxdy * pNode->dUdP; 566 *(pNode->fPN) += dxdy * pNode->dUdN; 567 } 568 } 569 } 570 571 /* Handle neighbor and edge dependent terms */ 572 pNode = pElem->pTLNode; 573 if ( pNode->nodeType != CONTACT ) { 574 *(pNode->fPsiPsiiP1) -= dyOverDx; 575 *(pNode->fPsiPsijP1) -= dxOverDy; 576 if ( pElem->elemType == SEMICON ) { 577 *(pNode->fNN) += dy * pTEdge->dJnDn + dx * pLEdge->dJnDn; 578 *(pNode->fPP) += dy * pTEdge->dJpDp + dx * pLEdge->dJpDp; 579 *(pNode->fNPsiiP1) += dy * pTEdge->dJnDpsiP1; 580 *(pNode->fNNiP1) += dy * pTEdge->dJnDnP1; 581 *(pNode->fPPsiiP1) += dy * pTEdge->dJpDpsiP1; 582 *(pNode->fPPiP1) += dy * pTEdge->dJpDpP1; 583 *(pNode->fNPsijP1) += dx * pLEdge->dJnDpsiP1; 584 *(pNode->fNNjP1) += dx * pLEdge->dJnDnP1; 585 *(pNode->fPPsijP1) += dx * pLEdge->dJpDpsiP1; 586 *(pNode->fPPjP1) += dx * pLEdge->dJpDpP1; 587 } 588 } 589 pNode = pElem->pTRNode; 590 if ( pNode->nodeType != CONTACT ) { 591 *(pNode->fPsiPsiiM1) -= dyOverDx; 592 *(pNode->fPsiPsijP1) -= dxOverDy; 593 if ( pElem->elemType == SEMICON ) { 594 *(pNode->fNN) += -dy * pTEdge->dJnDnP1 + dx * pREdge->dJnDn; 595 *(pNode->fPP) += -dy * pTEdge->dJpDpP1 + dx * pREdge->dJpDp; 596 *(pNode->fNPsiiM1) += dy * pTEdge->dJnDpsiP1; 597 *(pNode->fNNiM1) -= dy * pTEdge->dJnDn; 598 *(pNode->fPPsiiM1) += dy * pTEdge->dJpDpsiP1; 599 *(pNode->fPPiM1) -= dy * pTEdge->dJpDp; 600 *(pNode->fNPsijP1) += dx * pREdge->dJnDpsiP1; 601 *(pNode->fNNjP1) += dx * pREdge->dJnDnP1; 602 *(pNode->fPPsijP1) += dx * pREdge->dJpDpsiP1; 603 *(pNode->fPPjP1) += dx * pREdge->dJpDpP1; 604 } 605 } 606 pNode = pElem->pBRNode; 607 if ( pNode->nodeType != CONTACT ) { 608 *(pNode->fPsiPsiiM1) -= dyOverDx; 609 *(pNode->fPsiPsijM1) -= dxOverDy; 610 if ( pElem->elemType == SEMICON ) { 611 *(pNode->fNN) += -dy * pBEdge->dJnDnP1 - dx * pREdge->dJnDnP1; 612 *(pNode->fPP) += -dy * pBEdge->dJpDpP1 - dx * pREdge->dJpDpP1; 613 *(pNode->fNPsiiM1) += dy * pBEdge->dJnDpsiP1; 614 *(pNode->fNNiM1) -= dy * pBEdge->dJnDn; 615 *(pNode->fPPsiiM1) += dy * pBEdge->dJpDpsiP1; 616 *(pNode->fPPiM1) -= dy * pBEdge->dJpDp; 617 *(pNode->fNPsijM1) += dx * pREdge->dJnDpsiP1; 618 *(pNode->fNNjM1) -= dx * pREdge->dJnDn; 619 *(pNode->fPPsijM1) += dx * pREdge->dJpDpsiP1; 620 *(pNode->fPPjM1) -= dx * pREdge->dJpDp; 621 } 622 } 623 pNode = pElem->pBLNode; 624 if ( pNode->nodeType != CONTACT ) { 625 *(pNode->fPsiPsiiP1) -= dyOverDx; 626 *(pNode->fPsiPsijM1) -= dxOverDy; 627 if ( pElem->elemType == SEMICON ) { 628 *(pNode->fNN) += dy * pBEdge->dJnDn - dx * pLEdge->dJnDnP1; 629 *(pNode->fPP) += dy * pBEdge->dJpDp - dx * pLEdge->dJpDpP1; 630 *(pNode->fNPsiiP1) += dy * pBEdge->dJnDpsiP1; 631 *(pNode->fNNiP1) += dy * pBEdge->dJnDnP1; 632 *(pNode->fPPsiiP1) += dy * pBEdge->dJpDpsiP1; 633 *(pNode->fPPiP1) += dy * pBEdge->dJpDpP1; 634 *(pNode->fNPsijM1) += dx * pLEdge->dJnDpsiP1; 635 *(pNode->fNNjM1) -= dx * pLEdge->dJnDn; 636 *(pNode->fPPsijM1) += dx * pLEdge->dJpDpsiP1; 637 *(pNode->fPPjM1) -= dx * pLEdge->dJpDp; 638 } 639 } 640 } 641 642 /* Calculate the Inversion-Layer Mobility Dependent Terms in Jac. */ 643 if ( MobDeriv && SurfaceMobility ) { 644 for ( pCh = pDevice->pChannel; pCh != NULL; 645 pCh = pCh->next ) { 646 /* Find effective height of oxide element at interface. */ 647 if ( pCh->type%2 == 0 ) { /* Vertical slice */ 648 ds = pCh->pNElem->dy / pCh->pNElem->epsRel; 649 } else { /* Horizontal slice */ 650 ds = pCh->pNElem->dx / pCh->pNElem->epsRel; 651 } 652 pElem = pCh->pSeed; 653 nextIndex = (pCh->type + 2)%4; 654 while (pElem && pElem->channel == pCh->id) { 655 TWO_mobDeriv( pElem, pCh->type, ds ); 656 pElem = pElem->pElems[ nextIndex ]; 657 } 658 } /* endfor pCh != NULL */ 659 } /* endif MobDeriv and SurfaceMobility */ 660 } 661 662 /* load only the Rhs vector */ 663 void TWO_rhsLoad(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)664 TWO_rhsLoad(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info) 665 { 666 TWOelem *pElem; 667 TWOnode *pNode; 668 TWOedge *pHEdge, *pVEdge; 669 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 670 int index, eIndex; 671 double *pRhs = pDevice->rhs; 672 double dx, dy, dxdy, dyOverDx, dxOverDy; 673 double dPsiT, dPsiB, dPsiL, dPsiR; 674 double rhsN, rhsP; 675 double generation; 676 double nConc, pConc; 677 double perTime; 678 679 /* first compute the currents */ 680 TWO_commonTerms( pDevice, TRUE, tranAnalysis, info ); 681 682 /* find reciprocal timestep */ 683 if ( tranAnalysis ) { 684 perTime = info->intCoeff[0]; 685 } 686 687 /* zero the rhs vector */ 688 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) { 689 pRhs[ index ] = 0.0; 690 } 691 692 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 693 pElem = pDevice->elements[ eIndex ]; 694 695 dx = 0.5 * pElem->dx; 696 dy = 0.5 * pElem->dy; 697 dxdy = dx * dy; 698 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 699 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 700 701 pTEdge = pElem->pTopEdge; 702 pBEdge = pElem->pBotEdge; 703 pLEdge = pElem->pLeftEdge; 704 pREdge = pElem->pRightEdge; 705 dPsiT = pTEdge->dPsi; 706 dPsiB = pBEdge->dPsi; 707 dPsiL = pLEdge->dPsi; 708 dPsiR = pREdge->dPsi; 709 710 /* load for all i,j */ 711 for ( index = 0; index <= 3; index++ ) { 712 pNode = pElem->pNodes[ index ]; 713 if ( pNode->nodeType != CONTACT ) { 714 if ( index <= 1 ) { 715 pHEdge = pTEdge; 716 } else { 717 pHEdge = pBEdge; 718 } 719 if ( index == 0 || index == 3 ) { 720 pVEdge = pLEdge; 721 } else { 722 pVEdge = pREdge; 723 } 724 /* Add surface state charges. */ 725 pRhs[ pNode->psiEqn ] += dx * pHEdge->qf; 726 pRhs[ pNode->psiEqn ] += dy * pVEdge->qf; 727 if ( pElem->elemType == SEMICON ) { 728 nConc = pDevice->devState0 [pNode->nodeN]; 729 pConc = pDevice->devState0 [pNode->nodeP]; 730 pRhs[ pNode->psiEqn ] += dxdy * (pNode->netConc + pConc - nConc); 731 732 /* Handle generation terms */ 733 rhsN = - dxdy * pNode->uNet; 734 rhsP = dxdy * pNode->uNet; 735 if ( AvalancheGen ) { 736 generation = TWOavalanche( pElem, pNode ); 737 rhsN += dxdy * generation; 738 rhsP -= dxdy * generation; 739 } 740 pRhs[ pNode->nEqn ] -= rhsN; 741 pRhs[ pNode->pEqn ] -= rhsP; 742 743 /* Handle dXdT continuity terms */ 744 if ( tranAnalysis ) { 745 pRhs[ pNode->nEqn ] += dxdy * pNode->dNdT; 746 pRhs[ pNode->pEqn ] -= dxdy * pNode->dPdT; 747 } 748 } 749 } 750 } 751 752 /* Handle neighbor and edge dependent terms */ 753 pNode = pElem->pTLNode; 754 if ( pNode->nodeType != CONTACT ) { 755 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL; 756 if ( pElem->elemType == SEMICON ) { 757 pRhs[ pNode->nEqn ] -= dy * pTEdge->jn + dx * pLEdge->jn; 758 pRhs[ pNode->pEqn ] -= dy * pTEdge->jp + dx * pLEdge->jp; 759 } 760 } 761 pNode = pElem->pTRNode; 762 if ( pNode->nodeType != CONTACT ) { 763 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR; 764 if ( pElem->elemType == SEMICON ) { 765 pRhs[ pNode->nEqn ] -= -dy * pTEdge->jn + dx * pREdge->jn; 766 pRhs[ pNode->pEqn ] -= -dy * pTEdge->jp + dx * pREdge->jp; 767 } 768 } 769 pNode = pElem->pBRNode; 770 if ( pNode->nodeType != CONTACT ) { 771 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR; 772 if ( pElem->elemType == SEMICON ) { 773 pRhs[ pNode->nEqn ] -= -dy * pBEdge->jn - dx * pREdge->jn; 774 pRhs[ pNode->pEqn ] -= -dy * pBEdge->jp - dx * pREdge->jp; 775 } 776 } 777 pNode = pElem->pBLNode; 778 if ( pNode->nodeType != CONTACT ) { 779 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL; 780 if ( pElem->elemType == SEMICON ) { 781 pRhs[ pNode->nEqn ] -= dy * pBEdge->jn - dx * pLEdge->jn; 782 pRhs[ pNode->pEqn ] -= dy * pBEdge->jp - dx * pLEdge->jp; 783 } 784 } 785 } 786 } 787 788 /* 789 * computation of current densities, recombination rates, 790 * mobilities and their derivatives 791 */ 792 void TWO_commonTerms(TWOdevice * pDevice,BOOLEAN currentOnly,BOOLEAN tranAnalysis,TWOtranInfo * info)793 TWO_commonTerms(TWOdevice *pDevice, BOOLEAN currentOnly, 794 BOOLEAN tranAnalysis, TWOtranInfo *info) 795 { 796 TWOelem *pElem; 797 TWOedge *pEdge; 798 TWOnode *pNode; 799 int index, eIndex; 800 int nextIndex; /* index of node to find next element */ 801 double psi1, psi2, nC, pC, nP1, pP1; 802 double dPsiN, dPsiP; 803 double bPsiN, dbPsiN, bMPsiN, dbMPsiN; 804 double bPsiP, dbPsiP, bMPsiP, dbMPsiP; 805 double muN, dMuN, muP, dMuP, rDx, rDy; 806 double psi, nConc = 0.0, pConc = 0.0; 807 double cnAug, cpAug; 808 double eSurf = 0.0; /* For channel mobilities */ 809 double qInt = 0.0; 810 TWOchannel *pCh; 811 812 /* evaluate all node (including recombination) and edge quantities */ 813 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 814 pElem = pDevice->elements[ eIndex ]; 815 cnAug = pElem->matlInfo->cAug[ELEC]; 816 cpAug = pElem->matlInfo->cAug[HOLE]; 817 for ( index = 0; index <= 3; index++ ) { 818 if ( pElem->evalNodes[ index ] ) { 819 pNode = pElem->pNodes[ index ]; 820 if ( pNode->nodeType != CONTACT ) { 821 psi = pDevice->dcSolution[ pNode->psiEqn ]; 822 if ( pElem->elemType == SEMICON ) { 823 nConc = pDevice->dcSolution[ pNode->nEqn ]; 824 pConc = pDevice->dcSolution[ pNode->pEqn ]; 825 if ( Srh ) { 826 recomb(nConc, pConc, 827 pNode->tn, pNode->tp, cnAug, cpAug, pNode->nie, 828 &pNode->uNet, &pNode->dUdN, &pNode->dUdP); 829 } else { 830 pNode->uNet = 0.0; 831 pNode->dUdN = 0.0; 832 pNode->dUdP = 0.0; 833 } 834 } 835 } else { 836 /* a contact node */ 837 psi = pNode->psi; 838 if ( pElem->elemType == SEMICON ) { 839 nConc = pNode->nConc; 840 pConc = pNode->pConc; 841 } 842 } 843 844 /* store info in the state tables */ 845 pDevice->devState0 [pNode->nodePsi] = psi; 846 if ( pElem->elemType == SEMICON ) { 847 pDevice->devState0 [pNode->nodeN] = nConc; 848 pDevice->devState0 [pNode->nodeP] = pConc; 849 if ( tranAnalysis && (pNode->nodeType != CONTACT) ) { 850 pNode->dNdT = integrate( pDevice->devStates, info, pNode->nodeN ); 851 pNode->dPdT = integrate( pDevice->devStates, info, pNode->nodeP ); 852 } 853 } 854 } 855 } 856 for ( index = 0; index <= 3; index++ ) { 857 if ( pElem->evalEdges[ index ] ) { 858 pEdge = pElem->pEdges[ index ]; 859 pNode = pElem->pNodes[ index ]; 860 if ( pNode->nodeType != CONTACT ) { 861 psi1 = pDevice->dcSolution[pNode->psiEqn]; 862 } else { 863 psi1 = pNode->psi; 864 } 865 pNode = pElem->pNodes[ (index + 1) % 4 ]; 866 if ( pNode->nodeType != CONTACT ) { 867 psi2 = pDevice->dcSolution[pNode->psiEqn]; 868 } else { 869 psi2 = pNode->psi; 870 } 871 if ( index <= 1 ) { 872 pEdge->dPsi = psi2 - psi1; 873 } else { 874 pEdge->dPsi = psi1 - psi2; 875 } 876 pDevice->devState0 [pEdge->edgeDpsi] = pEdge->dPsi; 877 878 if ( pElem->elemType == SEMICON ) { 879 /* Calculate weighted driving forces - wdfn & wdfp for the edge */ 880 dPsiN = pEdge->dPsi + pEdge->dCBand; 881 dPsiP = pEdge->dPsi - pEdge->dVBand; 882 bernoulli( dPsiN, &bPsiN, &dbPsiN, 883 &bMPsiN, &dbMPsiN, !currentOnly ); 884 bernoulli( dPsiP, &bPsiP, &dbPsiP, 885 &bMPsiP, &dbMPsiP, !currentOnly); 886 if ( index <= 1 ) { 887 nC = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeN); 888 nP1 = *(pDevice->devState0 + pElem->pNodes[ index+1 ]->nodeN); 889 pC = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeP); 890 pP1 = *(pDevice->devState0 + pElem->pNodes[ index+1 ]->nodeP); 891 } else { 892 nC = *(pDevice->devState0 + pElem->pNodes[(index+1)%4]->nodeN); 893 nP1 = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeN); 894 pC = *(pDevice->devState0 + pElem->pNodes[(index+1)%4]->nodeP); 895 pP1 = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeP); 896 } 897 pEdge->wdfn = bPsiN * nP1 - bMPsiN * nC; 898 pEdge->wdfp = bPsiP * pC - bMPsiP * pP1; 899 pEdge->jn = 0.0; 900 pEdge->jp = 0.0; 901 if ( !currentOnly ) { 902 pEdge->dWnDpsiP1 = dbPsiN * nP1 - dbMPsiN * nC; 903 pEdge->dWnDn = - bMPsiN; 904 pEdge->dWnDnP1 = bPsiN; 905 pEdge->dWpDpsiP1 = dbPsiP * pC - dbMPsiP * pP1; 906 pEdge->dWpDp = bPsiP; 907 pEdge->dWpDpP1 = - bMPsiP; 908 pEdge->dJnDpsiP1 = 0.0; 909 pEdge->dJnDn = 0.0; 910 pEdge->dJnDnP1 = 0.0; 911 pEdge->dJpDpsiP1 = 0.0; 912 pEdge->dJpDp = 0.0; 913 pEdge->dJpDpP1 = 0.0; 914 } 915 } 916 } 917 } 918 } 919 920 /* DAG: calculate mobilities for channel elems */ 921 if ( SurfaceMobility ) { 922 for ( pCh = pDevice->pChannel; 923 pCh != NULL; pCh = pCh->next ) { 924 pElem = pCh->pNElem; 925 switch (pCh->type) { 926 case 0: 927 eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi ) 928 * pElem->epsRel / pElem->dy; 929 qInt = 0.5 * pElem->pBotEdge->qf; 930 break; 931 case 1: 932 eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi ) 933 * pElem->epsRel / pElem->dx; 934 qInt = 0.5 * pElem->pLeftEdge->qf; 935 break; 936 case 2: 937 eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi ) 938 * pElem->epsRel / pElem->dy; 939 qInt = 0.5 * pElem->pTopEdge->qf; 940 break; 941 case 3: 942 eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi ) 943 * pElem->epsRel / pElem->dx; 944 qInt = 0.5 * pElem->pRightEdge->qf; 945 break; 946 } 947 eSurf += qInt; 948 pElem = pCh->pSeed; 949 nextIndex = (pCh->type + 2)%4; 950 while (pElem && pElem->channel == pCh->id) { 951 TWO_mobility( pElem, eSurf ); 952 pElem = pElem->pElems[ nextIndex ]; 953 } 954 } /* endfor pCH != NULL */ 955 } /* endif SurfaceMobility */ 956 957 /* calculate the current densities assuming mobility value depend on Ex,Ey*/ 958 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 959 pElem = pDevice->elements[ eIndex ]; 960 961 rDx = 1.0 / pElem->dx; 962 rDy = 1.0 / pElem->dy; 963 for ( index = 0; index <= 3; index++ ) { 964 pEdge = pElem->pEdges[ index ]; 965 /* calculate conductive currents */ 966 if ( pElem->elemType == SEMICON ) { 967 /* get mobility for this edge */ 968 if ( !pElem->channel ) { 969 /* Calculate mobility for non-channel elements */ 970 muN = pElem->mun0; 971 dMuN = 0.0; 972 muP = pElem->mup0; 973 dMuP = 0.0; 974 dPsiN = pEdge->dPsi + pEdge->dCBand; 975 dPsiP = pEdge->dPsi - pEdge->dVBand; 976 if ( index%2 == 0 ) { 977 MOBfieldDep( pElem->matlInfo, ELEC, - dPsiN * rDx, &muN, &dMuN ); 978 MOBfieldDep( pElem->matlInfo, HOLE, - dPsiP * rDx, &muP, &dMuP ); 979 } else { 980 MOBfieldDep( pElem->matlInfo, ELEC, - dPsiN * rDy, &muN, &dMuN ); 981 MOBfieldDep( pElem->matlInfo, HOLE, - dPsiP * rDy, &muP, &dMuP ); 982 } 983 } else { 984 /* Retrieve previously calculated value. */ 985 muN = pElem->mun; 986 dMuN = 0.0; 987 muP = pElem->mup; 988 dMuP = 0.0; 989 } 990 switch ( index ) { 991 case 0: 992 muN *= pEdge->kPos * rDx; 993 dMuN *= pEdge->kPos * rDx * rDx; 994 muP *= pEdge->kPos * rDx; 995 dMuP *= pEdge->kPos * rDx * rDx; 996 break; 997 case 1: 998 muN *= pEdge->kNeg * rDy; 999 dMuN *= pEdge->kNeg * rDy * rDy; 1000 muP *= pEdge->kNeg * rDy; 1001 dMuP *= pEdge->kNeg * rDy * rDy; 1002 break; 1003 case 2: 1004 muN *= pEdge->kNeg * rDx; 1005 dMuN *= pEdge->kNeg * rDx * rDx; 1006 muP *= pEdge->kNeg * rDx; 1007 dMuP *= pEdge->kNeg * rDx * rDx; 1008 break; 1009 case 3: 1010 muN *= pEdge->kPos * rDy; 1011 dMuN *= pEdge->kPos * rDy * rDy; 1012 muP *= pEdge->kPos * rDy; 1013 dMuP *= pEdge->kPos * rDy * rDy; 1014 break; 1015 } 1016 /* Now that the mobility for this edge is known, do current */ 1017 pEdge->jn += muN * pEdge->wdfn; 1018 pEdge->jp += muP * pEdge->wdfp; 1019 if ( !currentOnly ) { 1020 pEdge->dJnDpsiP1 += muN * pEdge->dWnDpsiP1; 1021 pEdge->dJnDn += muN * pEdge->dWnDn; 1022 pEdge->dJnDnP1 += muN * pEdge->dWnDnP1; 1023 pEdge->dJpDpsiP1 += muP * pEdge->dWpDpsiP1; 1024 pEdge->dJpDp += muP * pEdge->dWpDp; 1025 pEdge->dJpDpP1 += muP * pEdge->dWpDpP1; 1026 if ( MobDeriv && (!pElem->channel) ) { 1027 pEdge->dJnDpsiP1 -= dMuN * pEdge->wdfn; 1028 pEdge->dJpDpsiP1 -= dMuP * pEdge->wdfp; 1029 } 1030 } 1031 } 1032 /* calculate displacement current only once */ 1033 if ( pElem->evalEdges[ index ] ) { 1034 if ( tranAnalysis ) { 1035 if ( index == 0 || index == 2 ) { 1036 /* horizontal edges */ 1037 pEdge->jd = -integrate(pDevice->devStates, info, 1038 pEdge->edgeDpsi) * rDx; 1039 } else { 1040 /* vertical edges */ 1041 pEdge->jd = -integrate(pDevice->devStates, info, 1042 pEdge->edgeDpsi) * rDy; 1043 } 1044 } 1045 } 1046 } 1047 } 1048 } 1049