Page MenuHomec4science
Files F65712463

gbar
No OneTemporary
Actions

File Metadata

Created
Wed, Jun 5, 17:18

gbar
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
nf = 0
first = 1
if output!= None:
f = open(output,'w')
for file in files:
nf = nf + 1
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)
# 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:
n = argmax(a)
if n==0 or n==len(r)-1:
# if the first is the max
xmx = r[n]
ymx = a[n]
#newpm = (p[0]+p[1]+p[2])/3.
### valeur moyenne de p ##
x = a*a*cos(p)
y = a*a*sin(p)
xx = sum(x[0:3])/sum(a[0:3]*a[0:3])
yy = sum(y[0:3])/sum(a[0:3]*a[0:3])
newpm = arctan2(yy,xx)
else:
x = r[n-1:n+2]
y = a[n-1:n+2]
A,xmx,ymx = fit_parabole(x,y)
### valeur moyenne de p ##
x = a*a*cos(p)
y = a*a*sin(p)
xx = sum(x[n-1:n+2])/sum(a[n-1:n+2]*a[n-1:n+2])
yy = sum(y[n-1:n+2])/sum(a[n-1:n+2]*a[n-1:n+2])
newpm = arctan2(yy,xx)
#newpm = (p[n-1]+p[n]+p[n+1])/3.
#newpm = -fmod(newpm/2.,2.*pi)
#newpm = newpm*180/pi
'''
# graph
g = SM.plot()
g.erase()
g.limits(r,a)
g.box()
g.connect(r,a)
g.show()
# graph
g = SM.plot()
g.erase()
g.limits(r,p)
g.box()
g.connect(r,p)
g.show()
'''
# graph
print r
print a
print p
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(newpm)
y = r*sin(newpm)
g.connect(x,y)
g.connect(-x,-y)
g.show()
if not first:
dpm = newpm - pm
if abs(dpm) > pi/2.:
dpm = dpm + pi
dt = newtm - tm
dpm = dpm/dt # conversion radian/Myr
# on regarde l'amplitude maximale, son rayon et l'angle
line = "%10.3f %10.5f %10.5f %10.5f %10.5f %02d" %(disk.tnow, xmx, ymx, dpm, newpm, npts)
print line
if output != None:
f.write(line+'\n')
else:
first = 0
pm = newpm
tm = newtm
if output != None:
f.close()

Event Timeline

c4science · Help