Skip to content
GitLab
Explore
Sign in
Adastra
meca_vol
Compare revisions
1c1d15bbcbe958e44641336266bad5cfe11f9079 to 3dc32fd1f1932a2d4b99a0e32c26163d1c6a6add
Commits on Source (2)
ajout de la fonction inv_ratioareasec_mach
· f079bc3e
pingu
authored
Sep 29, 2019
f079bc3e
ajout de la fonction inv_ratioareasec_mach
· 3dc32fd1
pingu
authored
Sep 29, 2019
3dc32fd1
Hide whitespace changes
Inline
Side-by-side
Engineering_tools.py
View file @
3dc32fd1
<<<<<<<
HEAD
import
numpy
as
np
from
numpy
import
linalg
as
LA
def
initial_vel
(
latitude
,
retrograde
=
False
):
"""
Initial velocity (m.s-1) depending on the launch point latitude (°)
"""
T_side
=
86164.10
# sidereal day (s)
R_T
=
6378e3
# earth raduis (m)
initial_velocity
=
(
2
*
np
.
pi
/
T_side
)
*
R_T
*
np
.
cos
(
np
.
deg2rad
(
np
.
abs
(
latitude
)))
if
retrograde
:
return
-
initial_velocity
else
:
return
initial_velocity
def
final_target_vel
(
h
,
a
):
"""
inputs:
h: injection altitude (m)
a: orbit semi major axis (m)
output:
orbit absolute injection velocity (m.s-1)
"""
mu
=
3.9860e14
# earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
R_T
=
6378e3
# earth raduis (m)
return
np
.
sqrt
(
2
*
mu
*
((
1
/
(
R_T
+
h
))
-
(
1
/
(
2
*
a
))))
def
staging_optim
(
isp1
,
isp2
,
eps1
,
eps2
,
Vc
,
Mp
):
"""
Optimise mass loading in the case of two stage rocket
inputs:
isp1, isp2 : specific impulse of the 1st and 2nd stage (s)
eps1, eps2 : structural coefficient of the 1st and 2nd stage (adimentional)
Vc : target critical velocity (m.s-1)
Mp : payload mass (Kg)
outputs:
M1 : mass of the 1st stage
M2 : mass of the 2nd stage
NB:
structure mass
structural coefficient = --------------------------------
structure mass + propellant mass
ctritical velocity = final velocity - initial velocity + losses
"""
# Notation :
#
# Mt = total mass
# q1 = M1/Mt
# q2 = M2/Mt
# qp = Mp/Mt
# omega = qp + q2
g0
=
9.81
# standard gravity (m.s-2)
omega_tab
=
np
.
linspace
(
0.001
,
0.9
,
100
)
# array of omega to test
qp_tab
=
[]
for
omega
in
omega_tab
:
q2
=
omega
+
eps1
*
(
1
-
omega
)
q2
=
np
.
power
(
q2
,
isp1
/
isp2
)
*
np
.
exp
(
Vc
/
(
g0
*
isp2
))
q2
=
1
/
q2
q2
=
1
-
q2
q2
*=
omega
/
(
1
-
eps2
)
qp
=
omega
-
q2
qp_tab
.
append
(
qp
)
qp
=
max
(
qp_tab
)
omega
=
omega_tab
[
qp_tab
.
index
(
qp
)]
q2
=
omega
-
qp
q1
=
1
-
q2
-
qp
M1
=
(
Mp
/
qp
)
*
q1
M2
=
(
Mp
/
qp
)
*
q2
return
M1
,
M2
def
standard_atmosphere
(
altitude
):
"""
Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model
input:
altitude (m)
outputs:
temperature (K)
pressure (Pa)
density (Kg.m-3)
"""
assert
altitude
<
84852
,
str
(
altitude
)
+
'
is higher than max altitude of the ISA model
'
g0
=
9.81
# standard gravity (m.s-2)
r
=
287
# specific gas constant (J.Kg-1.K-1)
# ISA model
h0_tab
=
np
.
array
([
0
,
11e3
,
20e3
,
32e3
,
47e3
,
51e3
,
71e3
,
84852
])
# base altitude (m)
a_tab
=
np
.
array
([
-
6.5e-3
,
0
,
1e-3
,
2.8e-3
,
0
,
-
2.8e-3
,
-
2.0e-3
])
# lapse rate (K.m-1)
# the base properties of each range are tabulated to reduce computational cost
T0_tab
=
np
.
array
([
288.15
,
216.65
,
216.65
,
228.65
,
270.65
,
270.65
,
214.65
])
# base temperature (K)
P0_tab
=
np
.
array
([
101325
,
22632
,
5474.9
,
868.02
,
110.91
,
66.939
,
3.9564
])
# base pressure (Pa)
index
=
0
while
altitude
>
h0_tab
[
index
+
1
]:
index
+=
1
a
=
a_tab
[
index
]
T0
=
T0_tab
[
index
]
P0
=
P0_tab
[
index
]
T
=
T0
+
a
*
(
altitude
-
h0_tab
[
index
])
# Toussaint equation
if
T
==
T0
:
# isohermal region
P
=
P0
*
np
.
exp
(
-
g0
*
(
altitude
-
h0_tab
[
index
])
/
(
r
*
T
))
else
:
# gradient region
P
=
P0
*
np
.
power
(
T
/
T0
,
-
g0
/
(
a
*
r
))
rho
=
P
/
(
r
*
T
)
return
T
,
P
,
rho
def
laval_nozzle
(
Pt
,
Tt
,
Pamb
,
Ath
,
Ae
,
gamma
,
r
,
tolerance
=
1e-5
):
"""
Compute the thrust assuming a started nozzle
inputs:
Pt : total pressure in the combustion chambre (Pa)
Tt : total temperature in the combustion chambre (K)
Pamb : ambient pressure (Pa)
Ath : throat section area (m2)
Ae : exit section area (m2)
gamma : ratio of specific heats (adimentional)
r : specific gas constant (J.Kg-1.K-1)
tolerance :
outputs:
Thrust (N)
Specific impulse (s)
N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area
"""
assert
gamma
>
1
,
'
gamma must be > 1
'
g0
=
9.81
# standard gravity (m.s-2)
A_ratio
=
lambda
M
,
gamma
:
(
1
/
M
)
*
np
.
power
((
2
/
(
gamma
+
1
))
*
(
1
+
((
gamma
-
1
)
/
2
)
*
M
*
M
),(
gamma
+
1
)
/
(
2
*
gamma
-
2
))
# compute the nozzle exit section area corresponding to the adapted regime
Madapted
=
np
.
sqrt
((
2
/
(
gamma
-
1
))
*
(
np
.
power
(
Pamb
/
Pt
,(
1
-
gamma
)
/
gamma
)
-
1
))
Aadapted
=
Ath
*
A_ratio
(
Madapted
,
gamma
)
if
Aadapted
>
Ae
:
# the nozzle is under-expanded
# compute the exit Mach number iteratively using dichotomy algorithm
A_ratio_target
=
Ae
/
Ath
Me_upper
=
Madapted
Me_lower
=
1
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
while
np
.
abs
(
error
)
>
tolerance
:
if
error
>
0
:
Me_lower
=
Me
else
:
Me_upper
=
Me
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
Te
=
Tt
/
(
1
+
((
gamma
-
1
)
/
2
)
*
Me
*
Me
)
# exit static temperature
Pe
=
Pt
*
np
.
power
(
Tt
/
Te
,
gamma
/
(
1
-
gamma
))
# exit static pressure
elif
Aadapted
<
Ae
:
# the nozzle is over-expanded
Mshock
=
1
step
=
0.1
A_ratio_target
=
Ae
/
Ath
Me_upper
=
1
Me_lower
=
0
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
while
np
.
abs
(
error
)
>
tolerance
:
if
error
<
0
:
Me_lower
=
Me
else
:
Me_upper
=
Me
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
Te
=
Tt
/
(
1
+
((
gamma
-
1
)
/
2
)
*
Me
*
Me
)
Pe
=
Pt
*
np
.
power
(
Tt
/
Te
,
gamma
/
(
1
-
gamma
))
while
(
Pe
-
Pamb
)
*
1e-5
>
tolerance
:
Mshock_old
=
Mshock
Mshock
+=
step
Pe_old
=
Pe
Pt_ratio
=
np
.
power
(((
gamma
+
1
)
*
Mshock
*
Mshock
)
/
((
gamma
-
1
)
*
Mshock
*
Mshock
+
2
),
gamma
/
(
gamma
-
1
))
*
np
.
power
((
gamma
+
1
)
/
(
2
*
gamma
*
Mshock
*
Mshock
-
gamma
+
1
),
1
/
(
gamma
-
1
))
Pi
=
Pt
*
Pt_ratio
Ai
=
Ath
/
Pt_ratio
A_ratio_target
=
Ae
/
Ai
Me_upper
=
1
Me_lower
=
0
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
while
np
.
abs
(
error
)
>
tolerance
:
if
error
<
0
:
Me_lower
=
Me
else
:
Me_upper
=
Me
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
Te
=
Tt
/
(
1
+
((
gamma
-
1
)
/
2
)
*
Me
*
Me
)
Pe
=
Pi
*
np
.
power
(
Tt
/
Te
,
gamma
/
(
1
-
gamma
))
if
(
Pe
>
Pamb
)
and
(
A_ratio
(
Mshock
,
gamma
)
>
Ae
/
Ath
):
# the normal shock is outside the divergent and all the flow inside is supersonic
A_ratio_target
=
Ae
/
Ath
Me_upper
=
Mshock
Me_lower
=
1
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
while
np
.
abs
(
error
)
>
tolerance
:
if
error
>
0
:
Me_lower
=
Me
else
:
Me_upper
=
Me
Me
=
(
Me_upper
+
Me_lower
)
/
2
error
=
A_ratio_target
-
A_ratio
(
Me
,
gamma
)
Te
=
Tt
/
(
1
+
((
gamma
-
1
)
/
2
)
*
Me
*
Me
)
# exit static temperature
Pe
=
Pt
*
np
.
power
(
Tt
/
Te
,
gamma
/
(
1
-
gamma
))
# exit static pressure
break
if
Pe
<
Pamb
:
Pe
=
Pe_old
Mshock
=
Mshock_old
step
/=
2
else
:
# the nozzle is adapted
Me
=
Madapted
Te
=
Tt
/
(
1
+
((
gamma
-
1
)
/
2
)
*
Me
*
Me
)
# exit static temperature
Pe
=
Pamb
Ve
=
Me
*
np
.
sqrt
(
gamma
*
r
*
Te
)
# exit velocity
# maximum mass flow rate in the case of a started nozzle
Q
=
np
.
power
(
2
/
(
gamma
+
1
),(
gamma
+
1
)
/
(
2
*
(
gamma
-
1
)))
*
np
.
sqrt
(
gamma
/
(
Tt
*
r
))
*
Pt
*
Ath
Thrust
=
Q
*
Ve
+
(
Pe
-
Pamb
)
*
Ae
isp
=
Thrust
/
(
Q
*
g0
)
return
Thrust
,
isp
def
adiabatic_flame
(
mech
,
T
,
P
,
X
):
"""
Compute adiatabic flame temperature
inputs:
mech : chemical mechanism in cti format
T : initial temperature (K)
P : pressure (Pa)
X : initial chemical composition
"""
# work in progress
def
injector
():
"""
"""
#%% Trajectory tools
def
solve_traj
():
"""
"""
# .......
#%% Astrodynamics tools
def
state_to_orbital_elem
(
r
,
v
):
"""
Convert state vector to orbital elements
inputs:
r : coordinate vector (m)
v : velocity vector (m.s-1)
output: {
'
a
'
,
'
e
'
,
'
i
'
,
'
capital_omega
'
,
'
small_omega
'
,
'
nu
'
}
where:
a : semi_major axis (m)
e : eccentricity (adimentional)
i : inclination (deg)
capital_omega : argument of the ascending node (deg)
small_omega : argument of the perigee (deg)
nu : true anomaly (deg)
NB: - An internal quantity (small) is considered for the case of undefined
classical orbital elements. When an element is undefined it is set to
zero and other(s) following it are measured according to this convention
- The function stops when r and v are aligned
- This function is widely inspired by the Matlab function written by
Prof. Josep Masdemart from UPC
"""
orbel
=
{
'
a
'
:
0
,
'
e
'
:
0
,
'
i
'
:
0
,
'
capital_omega
'
:
0
,
'
small_omega
'
:
0
,
'
nu
'
:
0
}
mu
=
3.9860e14
# earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
small
=
1e-10
h
=
np
.
cross
(
r
,
v
)
# angular momentum
rm
=
LA
.
norm
(
r
)
# module of the coordinate vector
vm
=
LA
.
norm
(
v
)
# module of the velocity vector
hm
=
LA
.
norm
(
h
)
# module of the angular momentum vector
if
hm
<
small
:
print
(
'
rectilinar or collision trajectory
'
)
return
n
=
np
.
array
([
-
h
[
1
],
h
[
0
],
0
])
nm
=
LA
.
norm
(
n
)
if
nm
<
small
:
n
=
np
.
array
([
1
,
0
,
0
])
nm
=
0
else
:
n
=
n
/
nm
nm
=
1
dotrv
=
np
.
dot
(
r
,
v
)
e
=
((
vm
*
vm
-
mu
/
rm
)
*
r
-
dotrv
*
v
)
/
mu
em
=
LA
.
norm
(
e
)
if
em
<
small
:
e
=
n
em
=
0
else
:
e
=
e
/
em
orbel
[
'
e
'
]
=
em
orbel
[
'
a
'
]
=
(
hm
*
hm
/
mu
)
/
(
1
-
em
*
em
)
aux
=
h
[
2
]
/
hm
if
np
.
abs
(
aux
)
>
1
:
aux
=
np
.
sign
(
aux
)
# all the angles are first computed in rad
orbel
[
'
i
'
]
=
np
.
arccos
(
aux
)
if
nm
>
small
:
aux
=
n
[
0
]
if
np
.
abs
(
aux
)
>
1
:
aux
=
np
.
sign
(
aux
)
orbel
[
'
capital_omega
'
]
=
np
.
arccos
(
aux
)
if
n
[
1
]
<
0
:
orbel
[
'
capital_omega
'
]
=
2
*
np
.
pi
-
orbel
[
'
capital_omega
'
]
if
em
>
small
:
aux
=
np
.
dot
(
e
,
n
)
if
np
.
abs
(
aux
)
>
1
:
aux
=
np
.
sign
(
aux
)
orbel
[
'
small_omega
'
]
=
np
.
arccos
(
aux
)
if
e
[
2
]
<
0
:
orbel
[
'
small_omega
'
]
=
2
*
np
.
pi
-
orbel
[
'
small_omega
'
]
aux
=
np
.
dot
(
e
,
r
)
/
rm
if
np
.
abs
(
aux
)
>
1
:
aux
=
np
.
sign
(
aux
)
orbel
[
'
nu
'
]
=
np
.
arccos
(
aux
)
if
em
>
small
:
if
dotrv
<
0
:
orbel
[
'
nu
'
]
=
2
*
np
.
pi
-
orbel
[
'
nu
'
]
else
:
ha
=
np
.
cross
(
e
,
r
)
if
np
.
dot
(
h
,
ha
)
<
0
:
orbel
[
'
nu
'
]
=
2
*
np
.
pi
-
orbel
[
'
nu
'
]
# convert the angles to deg
orbel
[
'
i
'
]
=
np
.
rad2deg
(
orbel
[
'
i
'
])
orbel
[
'
capital_omega
'
]
=
np
.
rad2deg
(
orbel
[
'
capital_omega
'
])
orbel
[
'
small_omega
'
]
=
np
.
rad2deg
(
orbel
[
'
small_omega
'
])
orbel
[
'
nu
'
]
=
np
.
rad2deg
(
orbel
[
'
nu
'
])
return
orbel
def
propagate_orbit
():
"""
"""
#%% Thermodynamic cycles
# the following functions provide the pressure, the mixuture composition
# and the initial temperature in the combustion chambre and properties of some
# components of the cycle.
def
pressure_fed
(
R
,
L
,
pi
,
Ltube
):
"""
input
p_fp: pressure of the fed pressure tank
p_co: pressure of the oxizider tank
p_ca: pressure of the fuel tank
d_fp: denisty of the fed pressure fluid
d_ox: density of the oxidizer
d_fu: density of the fuel
rajouter valves et ou et lesquelles ??
définir tank He
définir pipe He -> heat exchanger
définir pipe heat exchanger -> comb
définir pipe heat exchanger -> carb
définir tank comb
définir tank carb
définir pipe comb -> output
définir pipe carb -> output
"""
#vi = 0
t
=
0
dt
=
1e-1
tank_Hein
=
[
phe
,
vhe
,
0
,
dhe
,
rhe
,
0
,
rshe
,
lhe
,
che
]
tank_He
=
tank
(
tank_Hein
+
[
dt
])
tank_O2in
=
[
po
,
vo
,
do
,
co
,
ro
,
reo
]
tank_O2in
=
tank
(
p_fp
,
0
,
0
,
d_fp
,)
tank_CH4in
=
[
pch
,
vch
,
dch
,
cch
,
rch
,
rech
]
tank_CH4
=
tank
(
p_fp
,
0
,
0
,
d_fp
,)
pipeHeHe
=
pipes
(
p
,
v
,
d
,
mu
,
r
,
l
,
h
)
tank_
while
t
<
T
tank_He
=
tank
(
p_fp
,
0
,
0
,)
t
+=
dt
def
tank
(
p
,
v
,
d
,
dc
,
r
,
re
,
rs
,
l
,
c
,
dt
):
g
=
9.81
pi
=
3.1415926
"""
input
p : pression d
'
entrée
v : vitesse d
'
entrée
d : densité d
'
entrée
dc : densité du contenant
r : rayon
re : rayon d
'
entrée
rs : rayon de sortie
l : longueur
c : capacité de remplissage
dt : pas de temps
output
ps : pression sortie
vs : vitesse sortie
dc : densité
cs : capacité remplissage
"""
if
d_c
>=
0
pe
=
p
+
((
d
*
v
^
2
)
/
2
*
(
1
-
(
re
^
2
/
r
^
2
)))
ve
=
v
*
(
re
^
2
/
r
^
2
)
pm
=
pe
*
(
dc
/
d
)
+
((
dc
*
(
ve
^
2
)
/
2
)
*
(
1
-
(
d
/
dc
)
^
2
))
+
dc
*
l
*
g
*
(
1
-
2
*
c
)
vm
=
ve
*
(
d
/
dc
)
else
pm
=
p
-
(
d
*
g
*
h
)
-
(
8
*
mu
*
l
*
(
1
-
c
)
*
v
/
(
r
^
2
))
vm
=
v
ps
=
pm
+
((
dc
*
vm
^
2
)
/
2
*
(
1
-
(
r
^
2
/
rs
^
2
)))
vs
=
vm
*
(
r
^
2
/
rs
^
2
)
cs
=
c
-
((
vs
*
dc
*
dt
)
/
(
pi
*
r
^
2
*
l
))
return
ps
,
vs
,
dc
,
cs
def
pipes
(
p
,
v
,
d
,
mu
,
r
,
l
,
h
):
g
=
9.81
"""
input
p : pression d
'
entrée
v : vitesse d
'
entrée
d : densité d
'
entrée
mu : viscosité
r : rayon
l : longueur
h : longueur total projeté sur l
'
axe vertical
output
ps : pression sortie
vs : vitesse sortie
"""
pe
=
p
+
(
d
*
g
*
h
)
-
(
8
*
mu
*
l
*
v
/
(
r
^
2
))
ve
=
v
return
ps
,
vs
def
heat_exchanger
(
p
,
v
,
d
,
mu
,
r
,
l
,
L
):
g
=
9.81
"""
input
p : pression d
'
entrée
v : vitesse d
'
entrée
d : densité d
'
entrée
mu: viscosité dynamique
r : rayon
l : longueur
h : longueur total projeté sur l
'
axe vertical
L : chaleur latente de vaporisation
"""
pe
=
p
+
(
d
*
g
*
h
)
-
(
8
*
mu
*
l
*
v
/
(
r
^
2
))
ve
=
v
#heat_exchanger
ps
=
p
-
(
L
*
d
*
S
*
l
)
vs
=
ve
return
ps
,
vs
def
valves
(
x
,
p
,
v
,
dp
):
pi
=
3.1415926
"""
input
x: état imposé d
'
ouverture de la valve
p: pression d
'
entrée
v: vitesse d
'
entrée
dp: pressure drop
output
xs: état actuel d
'
ouverture de la valve
ps: pression de sortie de valves
vs: vitesse d
'
entrée
"""
ps
=
p
-
dp
#f(p,v,...)
vs
=
v
xs
=
x
return
xs
,
vs
,
ps
def
electric_turbo
():
"""
"""
def
closed_expander
():
"""
"""
def
dual_expander
():
"""
"""
def
open_expander
():
"""
"""
def
tap_off
():
"""
"""
def
gas_generator
():
"""
"""
def
ox_rich_staged_combust
():
"""
"""
def
full_flow_staged_combust
():
"""
"""
def
inv_ratioareasec_mach
(
rsa
,
gamma
,
subsup
)
"""
inputs
rsa : ratio des sections des aires
gamma : coefficient de Laplace
subsup : recherche d
'
un Mach sub ou supersonic (0|1)
Ce dernier input est lier à l
'
initialisation de la méthode.
"""
np
.
float128
gamma
=
1.4
epsilon
=
1e-10
if
subsup
=
1
:
x
=
5
else
:
x
=
0.2
iter
=
0
while
x
>
epsilon
and
iter
<
20
:
iter
+=
1
y
=
np
.
sqrt
((
2
/
(
gamma
+
1
)
*
(
1
+
(
gamma
-
1
)
*
(
x
**
2
)
/
2
))
**
((
gamma
+
1
)
/
(
gamma
-
1
))
/
x
**
2
)
-
rsa
xp
=
x
+
epsilon
yp
=
np
.
sqrt
((
2
/
(
gamma
+
1
)
*
(
1
+
(
gamma
-
1
)
*
(
xp
**
2
)
/
2
))
**
((
gamma
+
1
)
/
(
gamma
-
1
))
/
xp
**
2
)
-
rsa
x
=
x
-
epsilon
*
(
y
/
(
yp
-
y
))
#pour afficher chaque itération, décommenter
#print(x,y,iter)
return
x
=======
import
numpy
as
np
from
numpy
import
linalg
as
LA
...
...
@@ -564,3 +1126,33 @@ def ox_rich_staged_combust():
def
full_flow_staged_combust
():
"""
"""
def
inv_ratioareasec_mach
(
rsa
,
gamma
,
subsup
)
"""
inputs
rsa : ratio des sections des aires
gamma : coefficient de Laplace
subsup : recherche d
'
un Mach sub ou supersonic (0|1)
Ce dernier input est lier à l
'
initialisation de la méthode.
"""
np
.
float128
gamma
=
1.4
epsilon
=
1e-10
if
subsup
=
1
:
x
=
5
else
:
x
=
0.2
iter
=
0
while
x
>
epsilon
and
iter
<
20
:
iter
+=
1
y
=
np
.
sqrt
((
2
/
(
gamma
+
1
)
*
(
1
+
(
gamma
-
1
)
*
(
x
**
2
)
/
2
))
**
((
gamma
+
1
)
/
(
gamma
-
1
))
/
x
**
2
)
-
rsa
xp
=
x
+
epsilon
yp
=
np
.
sqrt
((
2
/
(
gamma
+
1
)
*
(
1
+
(
gamma
-
1
)
*
(
xp
**
2
)
/
2
))
**
((
gamma
+
1
)
/
(
gamma
-
1
))
/
xp
**
2
)
-
rsa
x
=
x
-
epsilon
*
(
y
/
(
yp
-
y
))
#pour afficher chaque itération, décommenter
#print(x,y,iter)
return
x
>>>>>>>
1
c1d15bbcbe958e44641336266bad5cfe11f9079