Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78481940
gbar3
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Aug 21, 03:52
Size
7 KB
Mime Type
text/x-python
Expires
Fri, Aug 23, 03:52 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20052710
Attached To
rGTOOLS Gtools
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
Log In to Comment