Page MenuHomec4science
Files F69463546

gbar3
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 1, 22:23

gbar3
View Options

This document is not UTF8. It was detected as Shift JIS and converted to UTF8 for display.
#!/usr/bin/env python
from Nbody import *
from Nbody import io
from numarray import *
from numarray import mlab as MLab
from numarray import fft as FFT
from Nbody.myNumeric import *
#from LinearAlgebra import *
#from numarray.LinearAlgebra import *
from numarray import linear_algebra as la
import sys
import copy
import SM
try:
from optparse import OptionParser
except ImportError:
from optik import OptionParser
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("-f",
action="store",
dest="output",
type="string",
default = None,
help="output file",
metavar=" FILE")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
####################################
def histogram(a,bins):
####################################
n = searchsorted(sort(a),bins)
n = concatenate([n,[len(a)]])
return n[1:]-n[:-1]
####################################
def tofrec(fmin,fmax,n,i):
####################################
nn = float(n)
if i<=n/2+1:
return (i/nn+0.5)*(fmax-fmin) + fmin
else:
return (i/nn-0.5)*(fmax-fmin) + fmin
####################################
def fourier(x,y):
####################################
'''
Fct = Sum_(m=1,n) amp_m cos( 2.pi f_m phi + phi_m )
= Sum_(m=1,n) amp_m cos( m phi + phi_m )
'''
dx = (mlab.max(x)-mlab.min(x))/len(x)
fmin = -1./(2.*dx)
fmax = -fmin
#f = fromfunction(lambda ii:tofrec(fmin,fmax,len(x),ii) ,(len(x),))
f = array([],Float)
for i in range(len(x)):
f = concatenate((f,tofrec(fmin,fmax,len(x),i)))
## FFT ##
ffft = fft.fft(y)
amp = sqrt(ffft.real**2 + ffft.imag**2) # amplitude du signal d'entree
phi = arctan2(ffft.imag,ffft.real) # phase a p/2 pres
phi = fmod(phi+2*pi,2.*pi) # de 0 a 2*pi (ok avec phase positive)
amp = 2.* amp / len(x) # doubler ok, /len(x) ???
f = f[0:len(x)/2] # on sort que les frequences positives
amp = amp[0:len(x)/2]
phi = phi[0:len(x)/2]
amp[0] = amp[0]/2.
return f,amp,phi
####################################
def parabole(a,x):
####################################
return a[0]*x**2 + a[1]*x + a[2]
####################################
def fit_parabole(x,y):
####################################
x12 = x[0]*x[0]
x22 = x[1]*x[1]
x32 = x[2]*x[2]
x1 = x[0]
x2 = x[1]
x3 = x[2]
A = array([[x12,x1,1],[x22,x2,1],[x32,x3,1]],Float)
B = y
a = la.solve_linear_equations(A,B)
xmx = -0.5*a[1]/a[0]
ymx = parabole(a,xmx)
return a,xmx,ymx
####################################
def khi2(r,p,a,angle):
####################################
'''
kind of khi2, weigted with the number of points
to a certain power
'''
e = 1/a
p = p-angle
x = r*cos(p)
y = r*sin(p)
# model y = 0
N = len(y)
#khi2 = sum( y**2/e**2 ) / N**6
#khi2 = sum( y**2/e**2 ) / (N-0.99)**5
khi2 = sum( y**2/e**2 ) / N**5
return khi2
####################################
def get_angle(r,p,a):
####################################
'''
algorithm used to find angle of the bar
'''
x = a*cos(p)
y = a*sin(p)
xx = sum(x)/sum(a)
yy = sum(y)/sum(a)
angle = arctan2(yy,xx)
return angle
####################################
def findbar(r,p,a):
####################################
'''
n is the number of point in the bar
'''
mn = 1e6
idx = 0
angle_opt = 0
n = 0
for i in range(3,len(r)+1):
rn = r[0:i]
pn = p[0:i]
an = a[0:i]
# estimation de l'angle
angle = get_angle(rn,pn,an)
# compute khi2
k = khi2(rn,pn,an,angle)
#print "%02d %02d %f"%(len(rn),len(r),k)
#print rn
if (k<mn):
idx = i
mn = k
angle_opt = angle
n = i
return angle_opt,n
######################################################################
# M A I N
######################################################################
# get options
files,options = parse_options()
ftype = options.ftype
output = options.output
if output!= None:
f = open(output,'w')
for file in files:
disk = Nbody(file,ftype=ftype)
nbody = disk.nbody
disk.histocenter()
# r : the radius used
# m : the modes computed
# m1 : the matrix of the amplitude
# m2 : the matrix of the phases
rs,m,m1,m2 = disk.dmodes(nr=32,nm=16,rm=16)
m1 = turnup(m1,0)
m2 = turnup(m2,0)
# plot
'''
i = 0
for r in rs:
g = SM.plot()
g.erase()
g.limits(0,10,m1[i,:])
g.box()
g.connect(m.astype(Float),m1[i,:])
g.show()
i=i+1
'''
r = []
a = []
p = []
m = []
inbar = 0
for i in range(len(rs)):
amp = m1[i,:] # extract amplitude
phi = m2[i,:] # extract phase
amx = argmax(amp[1:])+1
if not inbar and amx == 2:
inbar = 1
#print "in bar now"
if inbar and amx != 2:
inbar = 0
#print "out of bar now"
break
if inbar:
if phi[2]!=0: # avoid bad values (see dmodes)
r.append(rs[i])
a.append(amp[2])
p.append(phi[2])
m.append(amx)
###############################
#
m = array(m,Float)
newtm = disk.tnow
r = array(r,Float)
a = array(a,Float)
p = array(p,Float)
m = array(m,Float)
#print p
p = where((p<0),p+2*pi,p)
#print p*180/pi
#p = p/2. #!!! to get a correct spiral
#for i in range(len(r)):
# print r[i],p[i]
# recherche de la valeur maximale
### d騁ermination du maximum ###
npts = len(p)
if npts < 3:
print "found no bar at t=%8.3f (less than 3 pts)"%(disk.tnow)
else:
angle,n = findbar(r,p,a)
angle = angle/2 #!!! if not done before
rmx = max(r[0:n])
amx = max(a[0:n])
'''
x = r*cos(p)
y = r*sin(p)
g = SM.plot()
g.erase()
g.limits(-20,20,-20,20)
g.ptype(10,3)
g.box()
g.connect(x,y)
g.connect(-x,-y)
g.points(x,y)
g.ctype('red')
x = r[0:n]*cos(angle)
y = r[0:n]*sin(angle)
g.connect(x,y)
g.connect(-x,-y)
g.ctype('black')
g.show()
'''
line = "%10.3f %10.5f %10.5f %10.5f %02d" %(disk.tnow, angle, rmx, amx, n)
print line
if output != None:
f.write(line+'\n')
if output != None:
f.close()

Event Timeline

c4science · Help