Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85136212
gbar2
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
Fri, Sep 27, 00:49
Size
5 KB
Mime Type
text/x-python
Expires
Sun, Sep 29, 00:49 (2 d)
Engine
blob
Format
Raw Data
Handle
21133471
Attached To
rGTOOLS Gtools
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
Log In to Comment