octave_packages/secs2d-0.0.8/Utilities/Ubern.m

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