Page MenuHomec4science
Files F108733978

gbar2
No OneTemporary
Actions

File Metadata

Created
Fri, Apr 18, 05:21

gbar2
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
######################################################################
# 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 rr in r:
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)-1):
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:
r.append((rs[i]+rs[i+1])/2.) # radius
a.append(amp[amx]) # ampltude max mode
p.append(phi[amx]) # phase max mode
m.append(amx) # max mode
#print "add value"
###############################
#
m = array(m,Float)
newtm = disk.tnow
r = array(r,Float)
a = array(a,Float)
p = array(p,Float)
m = array(m,Float)
# calcul de l'angle de la barre :
npts = len(p)
if npts == 0:
print "found no bar at t=%8.3f (less than 3 pts)"%(disk.tnow)
else:
# moyenne pond駻馥
x = a*a*cos(p)
y = a*a*sin(p)
xx = sum(x)/sum(a*a)
yy = sum(y)/sum(a*a)
optp = arctan2(yy,xx)
n = argmax(a)
#optp = p[n]
rmx = r[n]
amx = a[n]
# graph
x = r*cos(p)
y = r*sin(p)
g = SM.plot()
g.erase()
g.limits(-15,15,-15,15)
g.ptype(10,3)
g.box()
g.connect(x,y)
g.connect(-x,-y)
g.points(x,y)
g.points(-x,-y)
x = r*cos(optp)
y = r*sin(optp)
g.connect(x,y)
g.connect(-x,-y)
g.show()
line = "%10.3f %10.5f %10.5f %10.5f %02d" %(disk.tnow, rmx, amx, optp, npts)
print line
if output != None:
f.write(line+'\n')
if output != None:
f.close()

Event Timeline

c4science · Help