octave_packages/secs1d-0.0.8/Utilities/Ubern.m

   1 function [bp,bn]=Ubern(x)
   2
   3 % [bp,bn]=Ubern(x)
   4 %
   5 % Bernoulli function
   6 % bp = B(x)=x/(exp(x)-1)
   7 % bn = B(-x)=x+B(x)
   8
   9
  10  ## This file is part of
  11   ##
  12   ## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
  13   ## -------------------------------------------------------------------
  14   ## Copyright (C) 2004-2007  Carlo de Falco
  15   ##
  16   ##
  17   ##
  18   ##  SECS1D is free software; you can redistribute it and/or modify
  19   ##  it under the terms of the GNU General Public License as published by
  20   ##  the Free Software Foundation; either version 2 of the License, or
  21   ##  (at your option) any later version.
  22   ##
  23   ##  SECS1D is distributed in the hope that it will be useful,
  24   ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
  25   ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  26   ##  GNU General Public License for more details.
  27   ##
  28   ##  You should have received a copy of the GNU General Public License
  29   ##  along with SECS1D; If not, see <http://www.gnu.org/licenses/>.
  30
  31
  32 xlim=1e-2;
  33 ax=abs(x);
  34
  35 %
  36 % Calcola la funz. di Bernoulli per x=0
  37 %
  38
  39 if (ax == 0)
  40    bp=1.;
  41    bn=1.;
  42    return
  43 end;
  44
  45 %
  46 % Calcola la funz. di Bernoulli per valori
  47 % asintotici dell'argomento
  48 %
  49
  50 if (ax > 80),
  51    if (x >0),
  52       bp=0.;
  53       bn=x;
  54       return
  55    else
  56       bp=-x;
  57       bn=0.;
  58       return
  59    end;
  60 end;
  61
  62 %
  63 % Calcola la funz. di Bernoulli per valori
  64 % intermedi dell'argomento
  65 %
  66
  67 if (ax > xlim),
  68    bp=x/(exp(x)-1);
  69    bn=x+bp;
  70    return
  71 else
  72
  73 %
  74 % Calcola la funz. di Bernoulli per valori
  75 % piccoli dell'argomento mediante sviluppo
  76 % di Taylor troncato dell'esponenziale
  77 %
  78    ii=1;
  79    fp=1.;
  80    fn=1.;
  81    df=1.;
  82    segno=1.;
  83    while (abs(df) > eps),
  84      ii=ii+1;
  85      segno=-segno;
  86      df=df*x/ii;
  87      fp=fp+df;
  88      fn=fn+segno*df;
  89      bp=1./fp;
  90      bn=1./fn;
  91    end;
  92    return
  93 end
  94
  95
  96 % Last Revision:
  97 % $Author: adb014 $
  98 % $Date: 2008-02-04 16:26:27 +0100 (man, 04 feb 2008) $