Aster/Cata/cataSTA9/Macro/gene_vari_alea_ops.py

   1 #@ MODIF gene_vari_alea_ops Macro  DATE 16/10/2007   AUTEUR REZETTE C.REZETTE
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2006  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20
  21 from math import sqrt,log,exp,pi,atan2,tan
  22
  23 def gene_vari_alea_ops(self,**args):
  24   self.set_icmd(1)
  25   return 0
  26
  27 def gene_vari_alea_init(self,d):
  28   from Utilitai.Utmess import  UTMESS
  29   a     =self.etape['BORNE_INF']
  30   moyen =self.etape['VALE_MOY' ]
  31   TYPE  =self.etape['TYPE']
  32   if self['INIT_ALEA']!=None :
  33      jump =self.etape['INIT_ALEA' ]
  34      self.iniran(jump)
  35   if   TYPE=='EXP_TRONQUEE' :
  36      b     =self.etape['BORNE_SUP']
  37      if (a>=b) :
  38          UTMESS('F','PROBA0_1',valr=[a,b])
  39      elif (moyen<=a)or(moyen>=b) :
  40          UTMESS('F','PROBA0_2',valr=[a,moyen,b])
  41      k=1./(moyen-a)
  42      if (exp(-b*k)<1.E-12) :
  43          UTMESS('F','PROBA0_3')
  44      # résolution par point fixe
  45      eps   =1.E-4
  46      nitmax=100000
  47      test  =0.
  48      while abs((test-k)/k)>eps :
  49          test = k
  50          k    = 1./(moyen-(a*exp(-a*k) - b*exp(-b*k))/(exp(-a*k) - exp(-b*k)))
  51      # génération de la variable aléatoire
  52      alpha = exp(-a*k) - exp(-b*k)
  53      self.sd.valeur=-( log(exp(-a*k)-alpha*self.getran()[0] ) ) /k
  54   elif TYPE=='EXPONENTIELLE' :
  55      if (moyen<=a) :
  56         UTMESS('F','PROBA0_4',valr=[moyen,a])
  57      v = moyen-a
  58      u=self.getran()[0]
  59      x = -log(1-u)
  60      self.sd.valeur=a + v*x
  61   elif TYPE=='GAMMA'         :
  62      delta =self.etape['COEF_VAR' ]
  63      if (moyen<=a) :
  64         UTMESS('F','PROBA0_4',valr=[moyen,a])
  65      v = moyen-a
  66      alpha = 1./(delta**2)
  67      if (alpha<=1.) :
  68         UTMESS('F','PROBA0_5')
  69      gamma2 = alpha-1.
  70      gamm1  = 1./gamma2
  71      beta   = sqrt(2.*alpha-1.)
  72      beta2  = 1./(beta**2)
  73      f0     = 0.5+(1./pi)*atan2(-gamma2/beta,1.)
  74      c1     = 1.-f0
  75      c2     = f0-0.5
  76      vref   = 0.
  77      vv     = -1.
  78 #
  79      while (-vv>vref) :
  80         u=self.getran()[0]
  81         gamdev = beta*tan(pi*(u*c1+c2))+gamma2
  82         unif=self.getran()[0]
  83         if unif<0. :
  84            UTMESS('F','PROBA0_6')
  85         vv= -log(unif)
  86         vref = log(1+beta2*((gamdev-gamma2)**2))+gamma2*log(gamdev*gamm1)-gamdev+gamma2
  87 #
  88      if vv<=0. :
  89         UTMESS('F','PROBA0_7')
  90      self.sd.valeur = a + v*(delta**2)*gamdev