]> Creatis software - CreaPhase.git/blob - octave_packages/secs1d-0.0.8/Utilities/Uscharfettergummel.m
Add a useful package (from Source forge) for octave
[CreaPhase.git] / octave_packages / secs1d-0.0.8 / Utilities / Uscharfettergummel.m
1 function A=Uscharfettergummel(nodes,Nnodes,elements,Nelements,acoeff,bcoeff,v)
2   
3 %A=Uscharfettergummel(nodes,Nnodes,elements,Nelements,acoeff,bcoeff,v)
4 %
5 % Builds the Scharfetter-Gummel  matrix for the 
6 % the discretization of the LHS 
7 % of the Drift-Diffusion equation:
8 %
9 % $ -(a(x) (u' - b v'(x) u))'= f $
10 %
11 % where a(x) is piecewise constant
12 % and v(x) is piecewise linear, so that 
13 % v'(x) is still piecewise constant
14 % b is a constant independent of x
15 % and u is the unknown
16 %
17
18   ## This file is part of 
19   ##
20   ## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
21   ## -------------------------------------------------------------------
22   ## Copyright (C) 2004-2007  Carlo de Falco
23   ##
24   ##
25   ##
26   ##  SECS2D is free software; you can redistribute it and/or modify
27   ##  it under the terms of the GNU General Public License as published by
28   ##  the Free Software Foundation; either version 2 of the License, or
29   ##  (at your option) any later version.
30   ##
31   ##  SECS2D is distributed in the hope that it will be useful,
32   ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
33   ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
34   ##  GNU General Public License for more details.
35   ##
36   ##  You should have received a copy of the GNU General Public License
37   ##  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
38
39   h=nodes(elements(:,2))-nodes(elements(:,1));
40   
41   c=acoeff./h;
42   Bneg=Ubernoulli(-(v(2:Nnodes)-v(1:Nnodes-1))*bcoeff,1);
43   Bpos=Ubernoulli( (v(2:Nnodes)-v(1:Nnodes-1))*bcoeff,1);
44   
45   
46   d0    = [c(1).*Bneg(1); c(1:end-1).*Bpos(1:end-1)+c(2:end).*Bneg(2:end); c(end)*Bpos(end)];
47   d1    = [1000;-c.* Bpos];
48   dm1   = [-c.* Bneg;1000];
49   
50   A = spdiags([dm1 d0 d1],-1:1,Nnodes,Nnodes);
51   
52   % Last Revision:
53   % $Author: adb014 $
54   % $Date: 2008-02-04 16:26:27 +0100 (man, 04 feb 2008) $