#!/usr/bin/python
# -*- coding: UTF-8 -*-
import numpy
from scipy import special
from cmath import *
from pylab import *
def sphjn(n,x):
return sqrt(pi/2./x)*special.jv(n+0.5,x)
def sphyn(n,x):
return sqrt(pi/2./x)*special.yv(n+0.5,x)
def sphh1n(n,x):
return sphjn(n,x)+1j*sphyn(n,x)
def sphh2n(n,x):
return sphjn(n,x)-1j*sphyn(n,x)
def MieCoeffSphPerfCond(a,ordre,freq):
c=3e8;
k0a=a*2./c*pi*freq;
A=numpy.zeros((freq.shape[0],ordre),complex)
B=numpy.zeros((freq.shape[0],ordre),complex)
for n in range(1,ordre+1):
N=float(2*n+1)/float(n)/float(n+1)
pj=pow(0-1j,n)
A[:,n-1] = -pj*N*sphjn(n,k0a)/sphh1n(n,k0a)
B[:,n-1] = (0-1j)*pj*N*(k0a*sphjn(n-1,k0a)-n*sphjn(n,k0a))/(k0a*sphh1n(n-1,k0a)-n*sphh1n(n,k0a))
return A, B
def sigBCSPerfCondSph(a,ordre,freq):
c=3e8
k0=2.*pi*freq/c
A=numpy.zeros((freq.shape[0],ordre),complex)
B=numpy.zeros((freq.shape[0],ordre),complex)
F0=numpy.zeros(freq.shape,complex)
A,B = MieCoeffSphPerfCond(a,ordre,freq)
for n in range(1,ordre+1):
pj=pow((0-1j),n-1)
F0=F0-pj*float(n)*float(n+1)*(A[:,n-1]+(0+1j)*B[:,n-1])/2.0
sigma=4*pi*abs(F0)*abs(F0)/k0/k0
return sigma
def test():
c=3e8
a=0.1
n=60
freq=numpy.arange(0.1e9,9.1e9,0.01e9)
k0a=a*2.0*pi*freq/c
sigma=sigBCSPerfCondSph(a,n,freq)/pi/a/a
plot(k0a, sigma, linewidth=1.0)
xlabel('k0a (Sphere circonference in wavelength)')
ylabel('BCS/pia2 Normalized echo area')
title("Backscattering Cross Section, metallic sphere")
grid(True)
show()
test()
SER d’une sphère parfaitement conductrice
13 décembre 2007 · Laisser un commentaire
Catégories : Python · Radars
Tagué : backscattering, BCS, conductor, scattering, SER, sphere
Différentiation automatique et Python
13 décembre 2007 · Laisser un commentaire
Voici un bout de code Python qui illustre le concept de différentiation automatique par surcharge d’opérateur. Le programme ci-dessous comporte des méthodes de test. On peut donc directement le donner à l’interpréteur. Pour plus d’informations.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import pylab
class adarray():
def __init__(self,v):
self.v = v
def __repr__(self):
return str(adarray(self.v))
def __str__(self):
return str(self.v)
def __add__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v+other.v)
def __radd__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v+other.v)
def __sub__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(self.v-other.v)
def __rsub__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
return adarray(other.v-self.v)
def __mul__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=other.v[:,0]*self.v[:,i]+self.v[:,0]*other.v[:,i]
return adarray(result)
def __rmul__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=other.v[:,0]*self.v[:,i]+self.v[:,0]*other.v[:,i]
return adarray(result)
def __div__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=self.v[:,0]/other.v[:,0]
o2=other.v[:,0]*other.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,0]*self.v[:,i]-self.v[:,0]*other.v[:,i])/o2
return adarray(result)
def __rdiv__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=other.v[:,0]/self.v[:,0]
o2=self.v[:,0]*self.v[:,0]
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,i]*self.v[:,0]-self.v[:,i]*other.v[:,0])/o2
return adarray(result)
def __pow__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=pow(self.v[:,0],other.v[:,0])
for i in range(1,self.v.shape[1]):
result[:,i]=(other.v[:,i]*log(self.v[:,i])+other.v[:,0]*self.v[:,i]/self.v[:,0])*result[:,0]
return adarray(result)
def __rpow__(self,other):
if not isinstance(other,adarray) :
o=numpy.zeros(self.v.shape)
o[:,0] = other
other = adarray(o)
result=numpy.zeros(self.v.shape)
result[:,0]=pow(other.v[:,0],self.v[:,0])
for i in range(1,self.v.shape[1]):
result[:,i]=(self.v[:,i]*log(other.v[:,i])+self.v[:,0]*other.v[:,i]/other.v[:,0])*result[:,0]
return adarray(result)
def sin(x):
if not isinstance(x,adarray) :
return numpy.sin(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.sin(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]*numpy.cos(x.v[:,0])
return adarray(result)
def cos(x):
if not isinstance(x,adarray) :
return numpy.sin(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.cos(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=-x.v[:,i]*numpy.sin(x.v[:,0])
return adarray(result)
def exp(x):
if not isinstance(x,adarray) :
return numpy.exp(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.exp(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]*numpy.exp(x.v[:,0])
return adarray(result)
def log(x):
if not isinstance(x,adarray) :
return numpy.log(x)
else :
result=numpy.zeros(x.v.shape)
result[:,0]=numpy.log(x.v[:,0])
for i in range(1,x.v.shape[1]):
result[:,i]=x.v[:,i]/x.v[:,0]
return adarray(result)
def CreateVar(base,calrange):
nbvar=calrange.shape[1]
for i in range(nbvar):
exec ("global %s%s" % (base,i))
exec ("globals()['%s%s'] = adarray(numpy.zeros([calrange.shape[0],calrange.shape[1]+1]))" % (base,i))
exec ("%s%s.v[:,0] = calrange[:,i]" % (base,i))
exec ("%s%s.v[:,i+1] = numpy.ones(calrange.shape[0])" % (base,i))
#exec ("print %s%s" % (base,i))
def testMonoDim():
### Cas monodimensionnel
plagex0=numpy.arange(0.001, 4., 0.1)
plage=numpy.array([plagex0]).T
CreateVar('x',plage)
## Opérations avec les scalaires
print "Opérations avec les scalaires :\n"
print "Matrice initiale : x = \n"+str(x0)
# Addition avec un scalaire après
print "Addition avec un scalaire après, \t x + 5 = \n%s" % (x0+5)
# Addition avec un scalaire avant
print "Addition avec un scalaire avant, \t 5 + x = \n%s" % (5+x0)
# Soustraction avec un scalaire après
print "Soustraction avec un scalaire après, \t x - 5 = \n%s" % (x0-5)
# Soustraction avec un scalaire avant
print "Soustraction avec un scalaire avant, \t 5 - x = \n%s" % (5-x0)
# Multiplication avec un scalaire après
print "Multiplication avec un scalaire après, \t x * 5 = \n%s" % (x0*5)
# Multiplication avec un scalaire avant
print "Multiplication avec un scalaire avant, \t 5 * x = \n%s" % (5*x0)
# Division avec un scalaire après
print "Division avec un scalaire après, \t x / 3 = \n%s" % (x0/3)
# Division avec un scalaire avant
print "Division avec un scalaire avant, \t 3 / x = \n%s" % (3/x0)
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t x^2 = \n%s" % (pow(x0,2))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t x^3 = \n%s" % (pow(x0,3))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t 2.^x = \n%s" % (pow(2,x0))
# Puissance avec un scalaire
print "Puissance avec un scalaire, \t 3.^x = \n%s" % (pow(3,x0))
## Opérations avec la variable
print "\nOpérations avec la variable :\n"
print "Vecteur initial : x = \n"+str(x0)
# Addition avec un pair
print "Addition,\t x + x = \n%s" % (x0+x0)
# Soustraction avec un pair
print "Soustraction,\t x - x = \n%s" % (x0-x0)
# Multiplication avec un pair
print "Multiplication,\t x * x = \n%s" % (x0*x0)
# Division avec un pair
print "Division,\t x / x = \n%s" % (x0/x0)
# Cosinus avec un pair
print "Cosinus,\t cos(x) = \n%s" % (cos(x0))
# Sinus avec un pair
print "Sinus, \t sin(x) = \n%s" % (sin(x0))
# Log avec un pair
print "Log, \t log(x) = \n%s" % (log(x0))
# Exp avec un pair
print "Exp, \t exp(x) = \n%s" % (exp(x0))
# Puissance avec un pair
print "Puissance avec un pair, \t x^x = \n%s" % (pow(x0,x0))
# Puissance avec un pair
print "Puissance avec un pair, \t x^(3*x) = \n%s" % (pow(x0,3*x0))
## Opérations complexes
print "\nOpérations complexes :\n"
print "Vecteur initial : x = \n"+str(x0)
# Addition
print "Addition,\t 10*x + 3*x = \n%s" % (10*x0+3*x0)
# Soustraction
print "Soustraction,\t 4*x - 5*x = \n%s" % (4*x0-5*x0)
# Multiplication
print "Multiplication,\t 3*x*4*x = \n%s" % (3*x0*4*x0)
# Division
print "Division,\t x/(x+1) = \n%s" % (x0/(x0+1))
# Cosinus
print "Cosinus,\t cos(pi*x/2-pi/2) = \n%s" % (cos(numpy.pi*x0/2-numpy.pi/2))
# Sinus
print "Sinus, \t sin(pi*x/2-pi/2) = \n%s" % (sin(numpy.pi*x0/2-numpy.pi/2))
# Log
print "Log, \t log(10*x)/log(10) = \n%s" % (log(10*x0)/log(10))
# Exp
print "Exp, \t x*exp(x/2) = \n%s" % (x0*exp(x0/2))
# Formule
print "Exp, \t x^5*exp(x/2)*sin(x0) = \n%s" % (pow(x0,5)*sin(x0)*exp(x0/2))
def testMonoDimGraph():
## Graphique
min=0.001
max=numpy.pi
nbpts=300
step=(max-min)/nbpts
print "Le pas vaut %s" %step
t=numpy.arange(min,max,step)
plage=numpy.array([t]).T
res=numpy.zeros(nbpts)
dres=numpy.zeros(nbpts)
CreateVar('x',plage)
y=sin(x0)*5*x0-1/(1+x0*x0)
pylab.figure(1)
pylab.plot(t,y.v[:,0])
pylab.plot(t,y.v[:,1])
pylab.show()
def testBiDim():
### Cas bidimensionnel
nbpts=6.
minx=1.
maxx=3.
rx=numpy.arange(minx, maxx, (maxx-minx)/nbpts)
miny=0.5
maxy=2.33
ry=numpy.arange(miny, maxy, (maxy-miny)/nbpts)
plage=numpy.array([rx,ry]).T
CreateVar('x',plage)
## Opérations avec les scalaires
print "Opérations avec les scalaires :\n"
print "Matrice initiale : x = \n"+str(x0)
print "Matrice initiale : y = \n"+str(x1)
# Addition avec un scalaire après
print "Addition avec un scalaire après, \t x + y + 5 = \n%s" % (x0+x1+5)
# Addition avec un scalaire avant
print "Addition avec un scalaire avant, \t 5 + x + y = \n%s" % (5+x0+x1)
# Soustraction avec un scalaire après
print "Soustraction avec un scalaire après, \t x + y - 5 = \n%s" % (x0+x1-5)
# Soustraction avec un scalaire avant
print "Soustraction avec un scalaire avant, \t 5 - x -y = \n%s" % (5-x0-x1)
# Multiplication avec un scalaire après
print "Multiplication avec un scalaire après, \t x * y * 5 = \n%s" % (x0*x1*5)
# Multiplication avec un scalaire avant
print "Multiplication avec un scalaire avant, \t 5 * x * y = \n%s" % (5*x0*x1)
# Division avec un scalaire après
print "Division avec un scalaire après, \t x * y / 3 = \n%s" % (x0*x1/3)
# Division avec un scalaire avant
print "Division avec un scalaire avant, \t 3 / (x * y) = \n%s" % (3/(x0*x1))
## Opérations complexes
print "\nOpérations complexes :\n"
print "Vecteur initial : x = \n"+str(x0)
print "Vecteur initial : y = \n"+str(x1)
# Addition
print "Addition,\t 10*x + 3*y = \n%s" % (10*x0+3*x1)
# Soustraction
print "Soustraction,\t 4*x - 5*y = \n%s" % (4*x0-5*x1)
# Multiplication
print "Multiplication,\t 3*x*4*y = \n%s" % (3*x0*4*x1)
# Division
print "Division,\t x/(y+1) = \n%s" % (x0/(x1+1))
# Cosinus
print "Cosinus,\t cos(pi*x*y/2) = \n%s" % (cos(numpy.pi*x0*x1/2))
# Sinus
print "Sinus, \t sin(pi*x*y/2) = \n%s" % (sin(numpy.pi*x0*x1/2))
# Log
print "Log, \t log(10*x*y)/log(10) = \n%s" % (log(10*x0*x1)/log(10))
# Exp
print "Exp, \t x*y*exp(x*y/2) = \n%s" % (x0*x1*exp(x0*x1/2))
# Complexe
print "Formule \t x*y*cos(x)-y*log(x*y) = \n%s" % (x0*x1*cos(x0)-x1*log(x0*x1))
testMonoDim()
testMonoDimGraph()
testBiDim()
Catégories : Python
Tagué : AD, automatic differentiation
Synchroniser une clé usb avec unison
13 décembre 2007 · Laisser un commentaire
Pour faire en sorte de ramener un système de fichiers identique au disque dur à la maison, on peut utiliser rsync ou unison. L’intérêt d’unison est qu’il fonctionne dans les deux sens et avec des systèmes microsoft et unix. La commande peut s’écrire :
unison -auto -batch ~/Documents /media/cleusb/backup
Les options auto et batch font qu’unison choisit ses propositions par défaut sans poser de questions interactives.
Faire un cron avec cette commande pour la lancer toutes les 5 minutes :
crontab -e
*/10 * * * 1-5 `unison -auto -batch ~/Documents /media/cleusb/backup`
