Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86920278
exercise_7.py
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, Oct 9, 10:24
Size
6 KB
Mime Type
text/x-python
Expires
Fri, Oct 11, 10:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21506270
Attached To
R9490 Homework_sp4e_Peruzzo_SáezUribe
exercise_7.py
View Options
#!/bin/env python3
import
main
import
numpy
as
np
import
argparse
from
scipy
import
optimize
def
readPositions
(
planet_name
,
directory
,
nb_steps
):
"""
This function reads a csv file with the state of the planet system and it create a matrix with the positions of a given planet at different times
:param planet_name: input the planet name with a NON capital first letter
:param directory: input the directory
:return: a numpy array of sizes nb_steps X 3
"""
m
=
nb_steps
n
=
3
planet_position
=
np
.
zeros
((
m
,
n
),
dtype
=
np
.
float64
)
# we need this if condition because the name of the files in the Trajectories folder is different from the name
# of the files that we get out of the code and are stored in the direcotry dumps
if
directory
==
"Trajectories/"
:
last
=
11
elif
directory
==
"dumps/"
:
last
=
10
else
:
print
(
"wrong directory for reference or solution files"
)
for
i
in
range
(
0
,
m
):
if
directory
==
"Trajectories/"
:
filename
=
directory
+
"step-"
+
'{:04}'
.
format
(
i
)
+
".csv"
elif
directory
==
"dumps/"
:
filename
=
directory
+
"step-"
+
'{:05}'
.
format
(
i
)
+
".csv"
else
:
print
(
"wrong directory for reference or solution files"
)
with
open
(
filename
,
'r'
)
as
f
:
data
=
f
.
readlines
()
for
line
in
data
:
words
=
line
.
split
()
if
words
[
last
]
==
planet_name
:
planet_position
[
i
,
0
:
n
]
=
np
.
array
(
words
[
0
:
n
],
dtype
=
np
.
float64
)
return
planet_position
def
computeError
(
positions_sol
,
positions_ref
):
"""
This function takes two matrices with the positions of a given planet at different times and it computes the integral error.
:param positions_sol: a numpy array of sizes 365 X 3 that contains the computed coordinates XYZ of the planets at different times
:param positions_ref: a numpy array of sizes 365 X 3 that contains the reference coordinates XYZ of the planets at different times
:return: type float64 containing the integral error of the simulation
"""
# todo: check if the numer of rows of positions_sol is equal to the one of positions_ref
error
=
np
.
sqrt
(
np
.
sum
(
np
.
power
(
np
.
subtract
(
positions_sol
,
positions_ref
),
2
)))
return
error
def
generateInput
(
scale
,
planet_name
,
input_filename
,
output_filename
):
"""
This function generates an input file from a given input file but by scaling velocity of a given planet.
:param scale: parameter to scale the velocity
:param planet_name: string with the name of the planet
:param input_filename: name of the file containing the input for the particle code to be scaled
:param output_filename: name of the file where to write the scaled input for the particle code
:return: file containing the input for the particle code for the scaled velocity
"""
out_file
=
open
(
output_filename
,
'w'
)
separator
=
' '
with
open
(
input_filename
,
'r'
)
as
f
:
data
=
f
.
readlines
()
for
line
in
data
:
words
=
line
.
split
()
if
words
[
10
]
==
planet_name
:
scale_velocity
=
np
.
multiply
(
scale
,
np
.
array
(
words
[
3
:
6
],
dtype
=
np
.
float64
))
string_velocity
=
np
.
array2string
(
scale_velocity
,
separator
=
' '
)
# There is a loss of precision here
out_file
.
write
(
separator
.
join
(
words
[
0
:
3
])
+
string_velocity
[
1
:
-
2
]
+
" "
+
separator
.
join
(
words
[
6
:
11
])
+
'
\n
'
)
else
:
out_file
.
write
(
separator
.
join
(
words
)
+
'
\n
'
)
out_file
.
close
()
return
out_file
def
launchParticles
(
input
,
nb_steps
,
freq
):
"""
This function launches the particle code for the specific case of planets and 1 day as time step, given an input
file, the number of steps and the frequency to dump
:param input: file containing the input for the particle code
:param nb_steps: number of steps
:param freq: frequency to dump
:return: nothing
"""
main
.
main
(
nb_steps
,
freq
,
input
.
name
,
"planet"
,
1
)
def
runAndComputeError
(
scale
,
planet_name
,
input
,
nb_steps
,
freq
):
"""
This function takes an input file containing the initial conditions of the planet system, and return the integral
error for a given scale of the velocity, a given planet, and a given number of steps
:param scale: parameter to scale the velocity
:param planet_name: string with the name of the planet
:param input: file containing the input for the particle code
:param nb_steps: number of steps
:param freq: frequency to dump
:return: type float64 containing the integral error of the simulation
"""
# creating the system data for the scaled velocity
input_scale
=
generateInput
(
scale
,
planet_name
,
input
.
name
,
"Init_scaled_onthefly.csv"
)
# launching the particle code for the scaled velocity
launchParticles
(
input_scale
,
nb_steps
,
freq
)
# getting the reference positions of the planet at different times
positions_ref
=
readPositions
(
planet_name
,
"Trajectories/"
,
nb_steps
)
# getting the computed positions of the planet at different times
positions_sol
=
readPositions
(
planet_name
,
"dumps/"
,
nb_steps
)
# computing the error
error
=
computeError
(
positions_sol
,
positions_ref
)
return
error
if
__name__
==
"__main__"
:
parser
=
argparse
.
ArgumentParser
(
description
=
'Compute_error'
)
parser
.
add_argument
(
'nsteps'
,
type
=
int
,
help
=
'specify the number of steps to perform'
)
parser
.
add_argument
(
'freq'
,
type
=
int
,
help
=
'specify the frequency for dumps'
)
parser
.
add_argument
(
'planet_name'
,
type
=
str
,
help
=
'specify the name of the planet'
)
parser
.
add_argument
(
'input_filename'
,
type
=
str
,
help
=
'input filename or the complete path to the file if it is not in the same folder as the script.'
)
args
=
parser
.
parse_args
()
# getting the parameters
planet_name
=
args
.
planet_name
scale
=
args
.
scale
nb_steps
=
args
.
nsteps
freq
=
args
.
freq
# getting the name of the file containing the input for the particle code to be scaled
input_filename
=
args
.
input_filename
# generate the input
input
=
open
(
input_filename
,
'r'
)
# compute the factor to scale the velocity of Mercury
minimum
=
optimize
.
fmin
(
runAndComputeError
,
1
,
args
=
(
planet_name
,
input
,
nb_steps
,
freq
))
# computing the error for a given scale
error
=
runAndComputeError
(
minimum
,
planet_name
,
input
,
nb_steps
,
freq
)
print
(
"The integral error between the reference and computed positions of the planet "
+
planet_name
+
" for a scale factor equal to "
,
scale
,
" is: "
)
print
(
error
)
Event Timeline
Log In to Comment