Aster/Cata/Utilitai/t_fonction.py

   1 #@ MODIF t_fonction Utilitai  DATE 31/05/2005   AUTEUR DURAND C.DURAND
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2005  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 from Numeric import *
  21 import copy
  22 import types
  23
  24 def interp(typ_i,val,x1,x2,y1,y2) :
  25   if typ_i==['LIN','LIN']: return y1+(y2-y1)*(val-x1)/(x2-x1)
  26   if typ_i==['LIN','LOG']: return exp(log(y1)+(val-x1)*(log(y2)-log(y1))/(x2-x1))
  27   if typ_i==['LOG','LOG']: return exp(log(y1)+(log(val)-log(x1))*(log(y2)-log(y1))/(log(x2)-log(x1)))
  28   if typ_i==['LOG','LIN']: return y1+(log(val)-log(x1))*(y2-y1)/(log(x2)-log(x1))
  29   if typ_i[0]=='NON'     :
  30                            if   val==x1 : return y1
  31                            elif val==x2 : return y2
  32                            else         : raise StandardError, 'fonction : interpolation NON'
  33 def is_ordo(liste) :
  34   listb=dict([(i,0) for i in liste]).keys()
  35   listb.sort()
  36   return liste==listb
  37
  38 class t_fonction :
  39   ### Classe pour fonctions réelles, équivalent au type aster = fonction_sdaster
  40   def __init__(self,vale_x,vale_y,para) :
  41     # création d'un objet fonction
  42     # vale_x et vale_y sont des listes de réels de meme longueur
  43     # para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
  44     pk=para.keys()
  45     pk.sort()
  46     if pk!=['INTERPOL','NOM_PARA','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
  47          raise StandardError, 'fonction : parametres incorrects'
  48     if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
  49          raise StandardError, 'fonction : parametre INTERPOL incorrect'
  50     if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
  51          raise StandardError, 'fonction : parametre PROL_DROITE incorrect'
  52     if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
  53          raise StandardError, 'fonction : parametre PROL_GAUCHE incorrect'
  54     self.vale_x    = array(vale_x)
  55     self.vale_y    = array(vale_y)
  56     self.para      = para
  57     if len(self.vale_x)!=len(self.vale_y) :
  58          raise StandardError, 'fonction : longueur abscisse <> longueur ordonnées'
  59     if not is_ordo(self.vale_x) :
  60          raise StandardError, 'fonction : abscisses non strictement croissantes'
  61
  62   def __add__(self,other) :
  63     # addition avec une autre fonction ou un nombre, par surcharge de l'opérateur +
  64     if   isinstance(other,t_fonction):
  65       para=copy.copy(self.para)
  66       vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
  67       fff=self.evalfonc(vale_x)
  68       ggg=other.evalfonc(vale_x)
  69       if   isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y+ggg.vale_y,para)
  70       else                              : return t_fonction(vale_x,fff.vale_y+ggg.vale_y,para)
  71     elif type(other) in [types.FloatType,types.IntType,types.ComplexType] :
  72       if   isinstance(self,t_fonction_c): return t_fonction_c(self.vale_x,self.vale_y+other,self.para)
  73       else                              : return t_fonction(self.vale_x,self.vale_y+other,self.para)
  74     else:  raise StandardError, 'fonctions : erreur de type dans __add__'
  75
  76   def __mul__(self,other) :
  77     # multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
  78     if   isinstance(other,t_fonction):
  79       para=copy.copy(self.para)
  80       vale_x,para['PROL_GAUCHE'],para['PROL_DROITE']=self.homo_support(other)
  81       fff=self.evalfonc(vale_x)
  82       ggg=other.evalfonc(vale_x)
  83       if   isinstance(self,t_fonction_c): return t_fonction_c(vale_x,fff.vale_y*ggg.vale_y,para)
  84       else                              : return t_fonction(vale_x,fff.vale_y*ggg.vale_y,para)
  85     elif type(other) in [types.FloatType,types.IntType] :
  86       return t_fonction(self.vale_x,self.vale_y*other,self.para)
  87     elif type(other) ==types.ComplexType :
  88       return t_fonction_c(self.vale_x,self.vale_y*other,self.para)
  89     else:  raise StandardError, 'fonctions : erreur de type dans __mul__'
  90
  91   def __repr__(self) :
  92     # affichage de la fonction en double colonne
  93     texte=[]
  94     for i in range(len(self.vale_x)) :
  95       texte.append('%f %f' % (self.vale_x[i],self.vale_y[i]))
  96     return '\n'.join(texte)
  97
  98   def __getitem__(self,other) :
  99     # composition de deux fonction F[G]=FoG=F(G(x))
 100     para=copy.copy(self.para)
 101     if other.para['NOM_RESU']!=self.para['NOM_PARA'] :
 102        raise StandardError,'''composition de fonctions : NOM_RESU1 et NOM_PARA2 incohérents '''
 103     para['NOM_PARA']==other.para['NOM_PARA']
 104     return t_fonction(other.vale_x,map(self,other.vale_y),para)
 105
 106   def __call__(self,val,tol=1.e-6):
 107     # méthode pour évaluer f(x)
 108     # tolérance, par défaut 1.e-6 en relatif sur la longueur de l'intervalle
 109     # adjacent, pour capter les erreurs d'arrondi en cas de prolongement exclu
 110     i=searchsorted(self.vale_x,val)
 111     n=len(self.vale_x)
 112     if i==0 :
 113       if self.para['PROL_GAUCHE']=='EXCLU'    :
 114          eps_g=(val-self.vale_x[0] )/(self.vale_x[1] -self.vale_x[0])
 115          if abs(eps_g)<=tol  : return self.vale_y[0]
 116          else                : raise StandardError, 'fonction évaluée hors du domaine de définition'
 117       else  :
 118          if self.para['PROL_GAUCHE']=='CONSTANT' : return self.vale_y[0]
 119          if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[0],
 120                                                                                            self.vale_x[1],
 121                                                                                            self.vale_y[0],
 122                                                                                            self.vale_y[1])
 123     elif i==n :
 124       if self.para['PROL_DROITE']=='EXCLU'    :
 125          eps_d=(val-self.vale_x[-1])/(self.vale_x[-1]-self.vale_x[-2])
 126          if abs(eps_d)<=tol  : return self.vale_y[-1]
 127          else                : raise StandardError, 'fonction évaluée hors du domaine de définition'
 128       else  :
 129          if self.para['PROL_DROITE']=='CONSTANT' : return self.vale_y[-1]
 130          if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val,self.vale_x[-1],
 131                                                                                            self.vale_x[-2],
 132                                                                                            self.vale_y[-1],
 133                                                                                            self.vale_y[-2])
 134     else :
 135       return interp(self.para['INTERPOL'],val,self.vale_x[i-1],
 136                                               self.vale_x[i],
 137                                               self.vale_y[i-1],
 138                                               self.vale_y[i])
 139
 140   def homo_support(self,other) :
 141     # renvoie le support d'abscisses homogénéisé entre self et other
 142     # i.e. si prolongement exclu, on retient plus grand min ou plus petit max, selon
 143     # si prolongement autorisé, on conserve les abscisses d'une fonction, extrapolantes
 144     # sur l'autre.
 145     # Pour les points intermédiaires : union et tri des valeurs des vale_x réunis.
 146     if other.vale_x[0]>self.vale_x[0]:
 147        if other.para['PROL_GAUCHE']!='EXCLU' : f_g=self
 148        else                                  : f_g=other
 149     else :
 150        if self.para['PROL_GAUCHE'] !='EXCLU' : f_g=other
 151        else                                  : f_g=self
 152     val_min    =f_g.vale_x[0]
 153     prol_gauche=f_g.para['PROL_GAUCHE']
 154     if self.vale_x[0]>other.vale_x[0]:
 155        if other.para['PROL_DROITE']!='EXCLU' : f_d=self
 156        else                                  : f_d=other
 157     else :
 158        if self.para['PROL_DROITE'] !='EXCLU' : f_d=other
 159        else                                  : f_d=self
 160     val_max    =f_d.vale_x[-1]
 161     prol_droite=f_d.para['PROL_DROITE']
 162     vale_x=concatenate((self.vale_x,other.vale_x))
 163     vale_x=clip(vale_x,val_min,val_max)
 164     vale_x=sort(dict([(i,0) for i in vale_x]).keys())
 165     return vale_x,prol_gauche,prol_droite
 166
 167   def cut(self,rinf,rsup,prec,crit='RELATIF') :
 168     # renvoie la fonction self dont on a 'coupé' les extrémités en x=rinf et x=rsup
 169     # pour la recherche de rinf et rsup dans la liste d'abscisses :
 170     #      prec=precision crit='absolu' ou 'relatif'
 171     para=copy.copy(self.para)
 172     para['PROL_GAUCHE']='EXCLU'
 173     para['PROL_DROITE']='EXCLU'
 174     if   crit=='ABSOLU' : rinf_tab=greater(abs(self.vale_x-rinf),prec)
 175     elif crit=='RELATIF': rinf_tab=greater(abs(self.vale_x-rinf),prec*rinf)
 176     else : raise StandardError, 'fonction : cut : critère absolu ou relatif'
 177     if   crit=='ABSOLU' : rsup_tab=greater(abs(self.vale_x-rsup),prec)
 178     elif crit=='RELATIF': rsup_tab=greater(abs(self.vale_x-rsup),prec*rsup)
 179     else : raise StandardError, 'fonction : cut : critère absolu ou relatif'
 180     if alltrue(rinf_tab) : i=searchsorted(self.vale_x,rinf)
 181     else                 : i=rinf_tab.tolist().index(0)+1
 182     if alltrue(rsup_tab) : j=searchsorted(self.vale_x,rsup)
 183     else                 : j=rsup_tab.tolist().index(0)
 184     vale_x=array([rinf,]+self.vale_x.tolist()[i:j]+[rsup,])
 185     vale_y=array([self(rinf),]+self.vale_y.tolist()[i:j]+[self(rsup),])
 186     return t_fonction(vale_x,vale_y,para)
 187
 188   def cat(self,other,surcharge) :
 189     # renvoie une fonction concaténée avec une autre, avec règles de surcharge
 190     para=copy.copy(self.para)
 191     if self.para['INTERPOL']!=other.para['INTERPOL'] : raise StandardError, 'concaténation de fonctions à interpolations différentes'
 192     if min(self.vale_x)<min(other.vale_x) :
 193             f1=self
 194             f2=other
 195     else                                  :
 196             f1=other
 197             f2=self
 198     para['PROL_GAUCHE']=f1.para['PROL_GAUCHE']
 199     if   surcharge=='GAUCHE' :
 200        i=searchsorted(f2.vale_x,f1.vale_x[-1])
 201        if i!=len(f2.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
 202        else                 : para['PROL_DROITE']=f1.para['PROL_DROITE']
 203        vale_x=array(f1.vale_x.tolist()+f2.vale_x[i:].tolist())
 204        vale_y=array(f1.vale_y.tolist()+f2.vale_y[i:].tolist())
 205     elif surcharge=='DROITE' :
 206        i=searchsorted(f1.vale_x,f2.vale_x[0])
 207        if i!=len(f1.vale_x) : para['PROL_DROITE']=f2.para['PROL_DROITE']
 208        else                 : para['PROL_DROITE']=f1.para['PROL_DROITE']
 209        vale_x=array(f1.vale_x[:i].tolist()+f2.vale_x.tolist())
 210        vale_y=array(f1.vale_y[:i].tolist()+f2.vale_y.tolist())
 211     return t_fonction(vale_x,vale_y,para)
 212
 213   def tabul(self) :
 214     # mise en forme de la fonction selon un vecteur unique (x1,y1,x2,y2,...)
 215     __tab=array([self.vale_x,self.vale_y])
 216     return ravel(transpose(__tab)).tolist()
 217
 218   def extreme(self) :
 219     # renvoie un dictionnaire des valeurs d'ordonnées min et max
 220     val_min=min(self.vale_y)
 221     val_max=max(self.vale_y)
 222     vm={}
 223     vm['min']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
 224                                                if self.vale_y[i]==val_min]
 225     vm['max']=[[self.vale_y[i],self.vale_x[i]] for i in range(len(self.vale_x))\
 226                                                if self.vale_y[i]==val_max]
 227     vm['min'].sort()
 228     vm['max'].sort()
 229     for item in vm['min'] : item.reverse()
 230     for item in vm['max'] : item.reverse()
 231     return vm
 232
 233   def trapeze(self,coef) :
 234     # renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
 235     trapz     = zeros(len(self.vale_y),Float)
 236     trapz[0]  = coef
 237     trapz[1:] = (self.vale_y[1:]+self.vale_y[:-1])/2*(self.vale_x[1:]-self.vale_x[:-1])
 238     prim_y=cumsum(trapz)
 239     para=copy.copy(self.para)
 240     para['PROL_GAUCHE']='EXCLU'
 241     para['PROL_DROITE']='EXCLU'
 242     if   para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
 243     elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
 244     return t_fonction(self.vale_x,prim_y,para)
 245
 246   def simpson(self,coef) :
 247     # renvoie la primitive de la fonction, calculée avec la constante d'intégration 'coef'
 248     para=copy.copy(self.para)
 249     para['PROL_GAUCHE']='EXCLU'
 250     para['PROL_DROITE']='EXCLU'
 251     if   para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='DEPL'
 252     elif para['NOM_RESU'][:4]=='ACCE' : para['NOM_RESU']='VITE'
 253     fm      = self.vale_y[0]
 254     fb      = self.vale_y[1]
 255     h2      = self.vale_x[1] - self.vale_x[0]
 256     tabl    = [coef,coef +(fb+fm)*h2/2.]
 257     prim_y  = copy.copy(tabl)
 258     iperm   = 0
 259     ip      = (1,0)
 260     eps     = 1.E-4
 261     for i in range(2,len(self.vale_x)) :
 262        h1  = h2
 263        h2  = self.vale_x[i] - self.vale_x[i-1]
 264        bma = h1 + h2
 265        fa  = fm
 266        fm  = fb
 267        fb  = self.vale_y[i]
 268        deltah = h2 - h1
 269        if h1==0. or h2==0. or abs( deltah / h1 ) <= eps :
 270           ct  = (1.,4.,1.)
 271        else :
 272           epsi = deltah / (h1*h2)
 273           ct   = (1.-epsi*h2,2.+epsi*deltah,1.+epsi*h1)
 274        tabl[iperm] = tabl[iperm] + (bma/6.)*(ct[0]*fa+ct[1]*fm+ct[2]*fb)
 275        prim_y.append(tabl[iperm])
 276        iperm       = ip[iperm]
 277     return t_fonction(self.vale_x,prim_y,para)
 278
 279   def derive(self) :
 280     # renvoie la dérivée de la fonction
 281     pas=self.vale_x[1:]-self.vale_x[:-1]
 282     pentes=(self.vale_y[1:]-self.vale_y[:-1])/(self.vale_x[1:]-self.vale_x[:-1])
 283     derive=(pentes[1:]*pas[1:]+pentes[:-1]*pas[:-1])/(pas[1:]+pas[:-1])
 284     derv_y=[pentes[0]]+derive.tolist()+[pentes[-1]]
 285     para=copy.copy(self.para)
 286     para['PROL_GAUCHE']='EXCLU'
 287     para['PROL_DROITE']='EXCLU'
 288     if   para['NOM_RESU'][:4]=='DEPL' : para['NOM_RESU']='VITE'
 289     elif para['NOM_RESU'][:4]=='VITE' : para['NOM_RESU']='ACCE'
 290     return t_fonction(self.vale_x,derv_y,para)
 291
 292   def inverse(self) :
 293     # renvoie l'inverse de la fonction
 294     # on intervertit vale_x et vale_y, on swape interpolation
 295     para=copy.copy(self.para)
 296     para['NOM_RESU']='TOUTRESU'
 297     para['NOM_PARA']=self.para['NOM_PARA']
 298     para['INTERPOL'].reverse()
 299     if para['PROL_GAUCHE']=='CONSTANT' : para['PROL_GAUCHE']='EXCLU'
 300     if para['PROL_DROITE']=='CONSTANT' : para['PROL_DROITE']='EXCLU'
 301     vale_x=self.vale_y
 302     vale_y=self.vale_x
 303     if not is_ordo(vale_x) :
 304        vale_x=vale_x[::-1]
 305        vale_y=vale_y[::-1]
 306     return t_fonction(vale_x,vale_y,para)
 307
 308   def abs(self) :
 309     # renvoie la mm fonction avec valeur absolue des ordonnées
 310     para=copy.copy(self.para)
 311     if para['PROL_GAUCHE']=='LINEAIRE' : para['PROL_GAUCHE']='EXCLU'
 312     if para['PROL_DROITE']=='LINEAIRE' : para['PROL_DROITE']='EXCLU'
 313     return t_fonction(self.vale_x,absolute(self.vale_y),para)
 314
 315   def evalfonc(self,liste_val) :
 316     # renvoie la mm fonction interpolée aux points définis par la liste 'liste_val'
 317     return t_fonction(liste_val,map(self,liste_val),self.para)
 318
 319   def sup(self,other) :
 320    # renvoie l'enveloppe supérieure de self et other
 321     para=copy.copy(self.para)
 322 #   commentaire : pour les prolongements et l'interpolation, c'est self
 323 #   qui prime sur other
 324     vale_x=self.vale_x.tolist()+other.vale_x.tolist()
 325 #   on ote les abscisses doublons
 326     vale_x=dict([(i,0) for i in vale_x]).keys()
 327     vale_x.sort()
 328     vale_x=array(vale_x)
 329     vale_y1=map(self, vale_x)
 330     vale_y2=map(other,vale_x)
 331     vale_y=map(max,vale_y1,vale_y2)
 332     return t_fonction(vale_x,vale_y,para)
 333
 334   def inf(self,other) :
 335    # renvoie l'enveloppe inférieure de self et other
 336     para=copy.copy(self.para)
 337 #   commentaire : pour les prolongements et l'interpolation, c'est self
 338 #   qui prime sur other
 339     vale_x=self.vale_x.tolist()+other.vale_x.tolist()
 340 #   on ote les abscisses doublons
 341     vale_x=dict([(i,0) for i in vale_x]).keys()
 342     vale_x.sort()
 343     vale_x=array(vale_x)
 344     vale_y1=map(self, vale_x)
 345     vale_y2=map(other,vale_x)
 346     vale_y=map(min,vale_y1,vale_y2)
 347     return t_fonction(vale_x,vale_y,para)
 348
 349   def suppr_tend(self) :
 350    # pour les corrections d'accélérogrammes
 351    # suppression de la tendance moyenne d'une fonction
 352     para=copy.copy(self.para)
 353     xy=sum(self.vale_x*self.vale_y)
 354     x0=sum(self.vale_x)
 355     y0=sum(self.vale_y)
 356     xx=sum(self.vale_x*self.vale_x)
 357     n=len(self.vale_x)
 358     a1 = ( n*xy - x0*y0) / (n*xx - x0*x0)
 359     a0 = (xx*x0 - x0*xy) / (n*xx - x0*x0)
 360     return t_fonction(self.vale_x,self.vale_y-a1*self.vale_x-a0,self.para)
 361
 362   def normel2(self) :
 363    # norme de la fonction
 364     __ex=self*self
 365     __ex=__ex.trapeze(0.)
 366     return sqrt(__ex.vale_y[-1])
 367
 368   def fft(self,methode) :
 369    # renvoie la transformée de Fourier rapide FFT
 370     import FFT
 371     para=copy.copy(self.para)
 372     para['NOM_PARA']='FREQ'
 373     if self.para['NOM_PARA']!='INST' :
 374        raise StandardError, 'fonction réelle : FFT : NOM_PARA=INST pour une transformée directe'
 375     pas = self.vale_x[1]-self.vale_x[0]
 376     for i in range(1,len(self.vale_x)) :
 377         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 378         if ecart>1.e-2 :
 379            raise StandardError, 'fonction réelle : FFT : la fonction doit etre à pas constant'
 380     n=int(log(len(self.vale_x))/log(2))
 381     if   methode=='TRONCATURE' :
 382        vale_y=self.vale_y[:2**n]
 383     elif methode=='PROL_ZERO'  :
 384        vale_y=self.vale_y.tolist()
 385        vale_y=vale_y+[0.]*(2**(n+1)-len(self.vale_x))
 386        vale_y=array(vale_y)
 387     vect=FFT.fft(vale_y)
 388     pasfreq =1./(pas*(len(vect)-1))
 389     vale_x  =[pasfreq*i for i in range(len(vect))]
 390     vale_y  =vect*pas
 391     return t_fonction_c(vale_x,vale_y,para)
 392
 393 class t_fonction_c(t_fonction) :
 394   ### Classe pour fonctions complexes, équivalent au type aster = fonction_c
 395   def tabul(self) :
 396     # mise en forme de la fonction selon un vecteur unique (x1,yr1,yi1,x2,yr2,yr2,...)
 397     __tab=array([self.vale_x,self.vale_y.real,self.vale_y.imag])
 398     return ravel(transpose(__tab)).tolist()
 399
 400   def __repr__(self) :
 401     # affichage de la fonction en double colonne
 402     texte=[]
 403     for i in range(len(self.vale_x)) :
 404       texte.append('%f %f + %f .j' % (self.vale_x[i],self.vale_y[i].real,self.vale_y[i].imag))
 405     return '\n'.join(texte)
 406
 407   def fft(self,methode,syme) :
 408    # renvoie la transformée de Fourier rapide FFT (sens inverse)
 409     import FFT
 410     para=copy.copy(self.para)
 411     para['NOM_PARA']='INST'
 412     if self.para['NOM_PARA']!='FREQ' :
 413        raise StandardError, 'fonction complexe : FFT : NOM_PARA=FREQ pour une transformée directe'
 414     pas = self.vale_x[1]-self.vale_x[0]
 415     for i in range(1,len(self.vale_x)) :
 416         ecart = abs(((self.vale_x[i]-self.vale_x[i-1])-pas)/pas)
 417         if ecart>1.e-3 :
 418            raise StandardError, 'fonction complexe : FFT : la fonction doit etre à pas constant'
 419     n=int(log(len(self.vale_x))/log(2))
 420     if   syme=='OUI' and len(self.vale_x)==2**n :
 421        vale_fonc=self.vale_y
 422     elif syme=='NON' and len(self.vale_x)==2**n :
 423        vale_fonc=self.vale_y.tolist()
 424        vale_fon1=vale_fonc[:]
 425        vale_fon1.reverse()
 426        vale_fonc=vale_fonc+vale_fon1
 427        vale_fonc=array(vale_fonc)
 428     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='PROL_ZERO' :
 429        vale_fonc=self.vale_y.tolist()+[complex(0.)]*(2**(n+1)-len(self.vale_x))
 430        vale_fon1=vale_fonc[:]
 431        vale_fon1.reverse()
 432        vale_fonc=vale_fonc+vale_fon1
 433        vale_fonc=array(vale_fonc)
 434     elif syme=='NON' and len(self.vale_x)!=2**n and methode=='TRONCATURE' :
 435        vale_fonc=self.vale_y[:2**n]
 436        vale_fonc=vale_fonc.tolist()
 437        vale_fon1=vale_fonc[:]
 438        vale_fon1.reverse()
 439        vale_fonc=vale_fonc+vale_fon1
 440        vale_fonc=array(vale_fonc)
 441     if   syme=='OUI' and len(self.vale_x)!=2**n :
 442        raise StandardError, 'fonction complexe : FFT : syme=OUI et nombre de points<>2**n'
 443     part1=vale_fonc[:len(vale_fonc)/2+1]
 444     part2=vale_fonc[1:len(vale_fonc)/2]
 445     part2=conjugate(part2)
 446     part2=part2.tolist()
 447     part2.reverse()
 448     vale_fonc=array(part1.tolist()+part2)
 449     vect=FFT.inverse_fft(vale_fonc)
 450     vect=vect.real
 451     pasfreq =1./(pas*(len(vect)-1))
 452     vale_x  =[pasfreq*i for i in range(len(vect))]
 453     pas2    =(1./self.vale_x[-1])*((len(self.vale_x))/float(len(vect)))
 454     vale_y  =vect/pas2
 455     return t_fonction(vale_x,vale_y,para)
 456
 457
 458 class t_nappe :
 459   ### Classe pour nappes, équivalent au type aster = nappe_sdaster
 460   def __init__(self,vale_para,l_fonc,para) :
 461     # création d'un objet nappe
 462     # vale_para est la liste de valeur des parametres (mot clé PARA dans DEFI_NAPPE)
 463     # para est un dictionnaire contenant les entrées PROL_DROITE, PROL_GAUCHE et INTERPOL (cf sd ASTER)
 464     # l_fonc est la liste des fonctions, de cardinal égal à celui de vale_para
 465     pk=para.keys()
 466     pk.sort()
 467     if pk!=['INTERPOL','NOM_PARA','NOM_PARA_FONC','NOM_RESU','PROL_DROITE','PROL_GAUCHE'] :
 468          raise StandardError, 'nappe : parametres incorrects'
 469     if para['INTERPOL'] not in [['NON','NON'],['LIN','LIN'],
 470                                 ['LIN','LOG'],['LOG','LOG'],['LOG','LIN'],] :
 471          raise StandardError, 'nappe : parametre INTERPOL incorrect'
 472     if para['PROL_DROITE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 473          raise StandardError, 'nappe : parametre PROL_DROITE incorrect'
 474     if para['PROL_GAUCHE'] not in ['EXCLU','CONSTANT','LINEAIRE'] :
 475          raise StandardError, 'nappe : parametre PROL_GAUCHE incorrect'
 476     self.vale_para    = array(vale_para)
 477     if type(l_fonc) not in (types.ListType,types.TupleType) :
 478          raise StandardError, 'nappe : la liste de fonctions fournie n est pas une liste'
 479     if len(l_fonc)!=len(vale_para) :
 480          raise StandardError, 'nappe : nombre de fonctions différent du nombre de valeurs du paramètre'
 481     for f in l_fonc :
 482       if not isinstance(f,t_fonction) and not isinstance(f,t_fonction_c) :
 483          raise StandardError, 'nappe : les fonctions fournies ne sont pas du bon type'
 484     self.l_fonc       = l_fonc
 485     self.para         = para
 486
 487   def __call__(self,val1,val2):
 488     # méthode pour évaluer nappe(val1,val2)
 489     i=searchsorted(self.vale_para,val1)
 490     n=len(self.vale_para)
 491     if i==0 :
 492       if val1==self.vale_para[0]  : return self.l_fonc[0](val2)
 493       if val1 <self.vale_para[0]  :
 494          if self.para['PROL_GAUCHE']=='EXCLU'    : raise StandardError, 'nappe évaluée hors du domaine de définition'
 495          if self.para['PROL_GAUCHE']=='CONSTANT' : return self.l_fonc[0](val2)
 496          if self.para['PROL_GAUCHE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 497                                                                  self.vale_para[0],
 498                                                                  self.vale_para[1],
 499                                                                  self.l_fonc[0](val2),
 500                                                                  self.l_fonc[1](val2))
 501     elif i==n :
 502       if val1==self.vale_para[-1] : return self.l_fonc[-1](val2)
 503       if val1 >self.vale_para[-1]  :
 504          if self.para['PROL_DROITE']=='EXCLU'    : raise StandardError, 'nappe évaluée hors du domaine de définition'
 505          if self.para['PROL_DROITE']=='CONSTANT' : return self.l_fonc[-1](val2)
 506          if self.para['PROL_DROITE']=='LINEAIRE' : return interp(self.para['INTERPOL'],val1,
 507                                                                  self.vale_para[-1],
 508                                                                  self.vale_para[-2],
 509                                                                  self.l_fonc[-1](val2),
 510                                                                  self.l_fonc[-2](val2))
 511     else :
 512       return interp(self.para['INTERPOL'],val1,self.vale_para[i-1],
 513                                                self.vale_para[i],
 514                                                self.l_fonc[i-1](val2),
 515                                                self.l_fonc[i](val2))
 516
 517   def __add__(self,other) :
 518     # addition avec une autre nappe ou un nombre, par surcharge de l'opérateur +
 519     l_fonc=[]
 520     if   isinstance(other,t_nappe):
 521       if self.vale_para!=other.vale_para : raise StandardError, 'nappes à valeurs de paramètres différentes'
 522       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other.l_fonc[i])
 523     elif type(other) in [types.FloatType,types.IntType] :
 524       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]+other)
 525     else:  raise StandardError, 't_nappe : erreur de type dans __add__'
 526     return t_nappe(self.vale_para,l_fonc,self.para)
 527
 528   def __mul__(self,other) :
 529     # multiplication avec une autre fonction ou un nombre, par surcharge de l'opérateur *
 530     l_fonc=[]
 531     if   isinstance(other,t_nappe):
 532       if self.vale_para!=other.vale_para : raise StandardError, 'nappes à valeurs de paramètres différentes'
 533       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other.l_fonc[i])
 534     elif type(other) in [types.FloatType,types.IntType] :
 535       for i in range(len(self.l_fonc)) : l_fonc.append(self.l_fonc[i]*other)
 536     else:  raise StandardError, 't_nappe : erreur de type dans __mul__'
 537     return t_nappe(self.vale_para,l_fonc,self.para)
 538
 539   def __repr__(self) :
 540     # affichage de la nappe en double colonne
 541     texte=[]
 542     for i in range(len(self.vale_para)) :
 543       texte.append('paramètre : %f' % self.vale_para[i])
 544       texte.append(repr(self.l_fonc[i]))
 545     return '\n'.join(texte)
 546
 547   def homo_support(self,other) :
 548     # renvoie la nappe self avec un support union de celui de self et de other
 549     # le support est la discrétisation vale_para et les discrétisations des fonctions
 550     if self==other : return self
 551     if self.para!=other.para : raise StandardError, 'combinaison de nappes à caractéristiques interpolation et prolongement différentes'
 552     vale_para=self.vale_para.tolist()+other.vale_para.tolist()
 553     vale_para=dict([(i,0) for i in vale_para]).keys()
 554     vale_para.sort()
 555     vale_para=array(vale_para)
 556     l_fonc=[]
 557     for val in vale_para :
 558         if   val in self.vale_para  : l_fonc.append(self.l_fonc[searchsorted(self.vale_para,val)])
 559         elif val in other.vale_para :
 560                                       other_fonc=other.l_fonc[searchsorted(other.vale_para,val)]
 561                                       new_vale_x=other_fonc.vale_x
 562                                       new_para  =other_fonc.para
 563                                       new_vale_y=[self(x) for x in new_vale_x]
 564                                       if isinstance(other_fonc,t_fonction)   :
 565                                                l_fonc.append(t_fonction(new_vale_x,new_vale_y,new_para))
 566                                       if isinstance(other_fonc,t_fonction_c) :
 567                                                l_fonc.append(t_fonction_c(new_vale_x,new_vale_y,new_para))
 568         else : raise StandardError, 'combinaison de nappes : incohérence'
 569     return t_nappe(vale_para,l_fonc,self.para)
 570
 571   def extreme(self) :
 572     # renvoie un dictionnaire des valeurs d'ordonnées min et max
 573     val_min=min([min(fonc.vale_y) for fonc in self.l_fonc])
 574     val_max=max([max(fonc.vale_y) for fonc in self.l_fonc])
 575     vm={'min':[],'max':[]}
 576     for i in range(len(self.vale_para)) :
 577         for j in range(len(self.l_fonc[i].vale_y)) :
 578           y = self.l_fonc[i].vale_y[j]
 579           if y==val_min : vm['min'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 580           if y==val_max : vm['max'].append([y,self.l_fonc[i].vale_x[j],self.vale_para[i]])
 581     vm['min'].sort()
 582     vm['max'].sort()
 583     for item in vm['min'] : item.reverse()
 584     for item in vm['max'] : item.reverse()
 585     return vm