Python planeta.Planeta類代碼示例

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt


condicion_inicial = [10., 0., 0., 0.3]

p = Planeta(condicion_inicial)

n_pasos=4000
dt=0.1

X=np.zeros(n_pasos)
Y=np.zeros(n_pasos)
VX=np.zeros(n_pasos)
VY=np.zeros(n_pasos)
ET=np.zeros(n_pasos)
t_values=np.zeros(n_pasos)

X[0]=condicion_inicial[0]
Y[0]=condicion_inicial[1]
VX[0]=condicion_inicial[2]
VY[0]=condicion_inicial[3]
ET[0]=p.energia_total()

p.avanza_rk4(dt)
X[1]=p.y_actual[0]
Y[1]=p.y_actual[1]

# -*- coding: utf-8 -*-

from planeta import Planeta
import matplotlib.pyplot as plt
import numpy as np

#Main

condicion_inicial = np.array([10.0, 0.0, 0.0, 0.25]) #Vy inicial = 0.25

p = Planeta(condicion_inicial)

#Inicalización
t_final = 900.0 #5 vueltas con euler
dt = 0.1
N_pasos= int(t_final/dt)

x = np.zeros(N_pasos)
y = np.zeros(N_pasos)
vx = np.zeros(N_pasos)
vy = np.zeros(N_pasos)

E = np.zeros(N_pasos)

[x[0],y[0],vx[0],vy[0]] = condicion_inicial
E[0] = p.energia_total()

#Integración
for i in range(1,N_pasos):
    p.avanza_euler(dt)
    xi, yi, vxi, vyi = p.y_actual

En este script se realizara el grafico de energia y trayectoria usando euler.
MIENTRAS MAS CHICO EL DT MAS PRECISO ES PARA LA ENERGIA PERO LA TRAYECTORIA NO

'''


'''
 Condiciones iniciales mas la creacion del objeto de clase Planeta.
 Paso arbitrario.
 Arreglos de ceros a los que se les ira anexando valores dado el metodo de euler.
'''

vy0=0.3
condicion_inicial = sp.array([10, 0, 0, vy0])
p = Planeta(condicion_inicial)
t_final =  6000
pasos = 50000
dt= t_final / (float)(pasos)
x = sp.zeros(pasos)
y = sp.zeros(pasos)
vx = sp.zeros(pasos)
vy = sp.zeros(pasos)
energia=sp.zeros(pasos)

[x[0],y[0],vx[0],vy[0]] = condicion_inicial
energia[0]=p.energia_total()

'''Caso euler explicito con alpha=0 .Se recorrera con un for para ir obteniendo para cada tiempo un vector de 4 componentes [x,y,vx,vy] '''

# -*- coding: utf-8 -*-


from planeta import Planeta
import matplotlib.pyplot as plt
import numpy as np

#Main

condicion_inicial = np.array([10.0, 0.0, 0.0, 0.25]) #Vy inicial = 0.25

p = Planeta(condicion_inicial)

#Inicalización
t_final = 900.0
dt = 0.1
N_pasos= int(t_final/dt)

x = np.zeros(N_pasos)
y = np.zeros(N_pasos)
vx = np.zeros(N_pasos)
vy = np.zeros(N_pasos)

E = np.zeros(N_pasos)

#Integración
[x[0],y[0],vx[0],vy[0]] = condicion_inicial
E[0] = p.energia_total()
for i in range(1,N_pasos):
    p.avanza_rk4(dt)
    xi, yi, vxi, vyi = p.y_actual

'''
Este script utiliza la clase planeta para integrar la trayectoria de
aproximadamente 5 orbitas con el metodo de verlet, graficando la trayectoria
y la energia vs tiempo en cada momento. Grafica dos veces energia vs tiempo
con el fin de dejarlo en dos escalas: una de cerca para ver a forma de la
funcion y otra mas de lejos para poder comparar con otros metodos menos
eficientes. Utiliza alpha=0 en el potencial.
'''
from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt

condicion_inicial = [10, 0, 0, 0.4]
p = Planeta(condicion_inicial)

N_steps = 8000
dt=4000./N_steps

t=np.linspace(0,4000,N_steps)

#creamos los arrays en que se guardara la informacion
x= np.zeros(N_steps)
y= np.zeros(N_steps)
vx= np.zeros(N_steps)
vy= np.zeros(N_steps)
energia = np.zeros(N_steps)

x[0]= 10
y[0]= 0
energia[0]=p.energia_actual

def orbitar(CI,solucion, prt=True):
    '''
    Esta función  es cool, porque le entregas las condiciones iniciales (CI)
    en forma de x0,y0,vx0,vy0, además de un string con el tipo de solución
    que quieres, como 'euler', 'rk4' o 'verlet' y nos devuelve las coordenadas
    x,y y la energia total del planeta. Además de un ploteo de la órbita, apropi-
    adamente labelado y un ploteo de la energía vs el tiempo.
    La ultima parte es para que plotee, el valor default es que
    plotee (true or false).
    Ejemplo de uso:
        orbitar([10,0,0,0.1],'verlet')
        x=[0,0.1,...]
        y=[0,0.3,...]
        energia=[7,8,...]
    '''
    if solucion=='euler':
        n= 2500 #grande y dsps se arregla
    if solucion=='rk4':
        n= 5000 #grande y dsps se arregla
    if solucion=='verlet':
        n=800
    Aiur= Planeta(CI) #se crea el planeta
    if solucion=='verlet_reloaded': #caso especial para la ultima parte
        n=6000
        Aiur= Planeta(CI,alpha=alpha2) #si es el verlet bkn se parcha
    dt= 1 #error si es muy chico
    tmax=n*dt
    x=[] #listas vacias a llenar
    y=[]
    pylon=[] #energia

    x= np.append(x,CI[0]) #llenamos con el primer termino
    y= np.append(y,CI[1])
    pylon= np.append(pylon,Aiur.energia_total())

    for i in range(tmax):   #aca se hace la iteracion del paso
        if solucion=='euler':
            Aiur.avanza_euler(dt)
        if solucion=='rk4':
            Aiur.avanza_rk4(dt)
        if solucion=='verlet' or solucion=='verlet_reloaded':
            Aiur.avanza_verlet(dt)
        xn,yn,vxn,vyn=Aiur.y_now
        x= np.append(x,xn)
        y= np.append(y,yn)
        pylon= np.append(pylon,Aiur.energia_total())
        times=np.linspace(0,tmax,tmax+1)
    if prt==True:
        fig1=plt.figure(1)  #aca se plotea, muy standar
        fig1.clf()
        if solucion=='euler':
            plt.plot(x,y,'r',label= 'Trayectoria metodo Euler')
            plt.title('Trayectoria metodo Euler')
        if solucion=='rk4':
            plt.plot(x,y,'r',label= 'Trayectoria metodo RK4')
            plt.title('Trayectoria metodo RK4')
        if solucion=='verlet' or solucion=='verlet_reloaded':
            plt.plot(x,y,'r',label= 'Trayectoria metodo Verlet')
            plt.title('Trayectoria metodo Verlet')

        plt.xlabel('X')
        plt.ylabel('Y')
        plt.savefig('verlet_rld.png')
        plt.draw()
        plt.show()

        fig2=plt.figure(2)
        fig2.clf()
        plt.plot(times, pylon,'g')
        plt.xlabel('Tiempo')
        plt.ylabel('Energia')
        plt.title('Tiempo vs Energia')
        plt.savefig('verlet_rld_energia.png')
        plt.draw()
        plt.show()
    return x,y,pylon

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt
dty0=0.3#condicion inicial de velocidad en y
condicion_inicial = np.array([10., 0, 0, dty0])

p = Planeta(condicion_inicial)

t_fin =  3000.#tiempo maximo
num_pasos = 2000+1#intervalos
dt= t_fin / (float)(num_pasos)#discretisacion
#crear los arreglos a usar mas adelante
x = np.zeros(num_pasos)
y = np.zeros(num_pasos)
dtx = np.zeros(num_pasos)
dty = np.zeros(num_pasos)
energia = np.zeros(num_pasos)
#se inician los valores para las variables
[x[0],y[0],dtx[0],dty[0]] = condicion_inicial
energia[0] = p.energia_total()

"""implementacion de verlet con alfa distinto de cero"""


"""
#da el primer paso usando RK 4, para continuar con el metodo de verlet
p.avanza_rk4(dt)
aux = p.y_actual
x[1] = aux[0]
y[1] = aux[1]
dtx[1] = aux[2]

# -*- coding: utf-8 -*-

'''
Programa que plotea las orbitas de un planeta. El metodo usado es verlet
'''

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt

fig=plt.figure()
fig.clf()

condicion_inicial = [10, 0, 0, 0.25]

Shakuras = Planeta(condicion_inicial)
N=3*np.int(1e4)
dt=1000./N
t=np.linspace(0,1000,N)
x=np.zeros(N)
y=np.zeros(N)
E=np.zeros(N)

x[0]=condicion_inicial[0]
y[0]=condicion_inicial[1]
E[0]=Shakuras.energia_total()

for n in range(1,N):
    Shakuras.avanza_verlet(dt)
    x[n]=Shakuras.y_actual[0]
    y[n]=Shakuras.y_actual[1]

#!/usr/bin/env python
# -*- coding: utf-8 -*-

'''
Script que utiliza Verlet para calcular y graficar aproximadamente
30 orbitas considerando alpha distinto de 0. Además grafica energia vs tiempo.
'''

import numpy as np
import matplotlib.pyplot as plt
from planeta import Planeta

condicion_inicial  =  [10, 0, 0, 0.3]
p  =  Planeta(condicion_inicial, alpha = 10**(-2.257))

# periodo aproximado por tercera ley de kepler
T = np.sqrt(4*np.pi**2*8**3) #  =  142 aproximadamente
t_total = 30*T
# se crean arreglos vacios de x, y, Energia
x = np.array([])
y = np.array([])
Energia = np.array([])
#Se agregan condiciones iniciales
x = np.append(x,10)
y = np.append(y,0)
Energia = np.append(Energia,p.energia_total())

dt = 1
t = 1
while t<= t_total:
    p.avanza_verlet(dt)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt


condicion_inicial = [10, 0, 0, 0.4]

p = Planeta(condicion_inicial)

dt = 0.1
iteraciones = 50000
x = np.zeros(iteraciones)
y = np.zeros(iteraciones)
energia = np.zeros(iteraciones)
tiempo = np.zeros(iteraciones)

for i in range(iteraciones):
    p.avanza_euler(dt)
    x[i] = p.y_actual[0]
    y[i] = p.y_actual[1]
    energia[i] =p.energia_total()
    tiempo [i] = p.t_actual

plt.figure(1)
plt.clf()
plt.plot(x,y,color='red')
plt.xlabel('x[m]')
plt.ylabel('y[m]')

from planeta import Planeta
import scipy.stats
import numpy as np


'''
PARTE3
Ahora considero el caso alfa no cero, mi rut es 189563752 y se usara verletvelocity para integrar la ecuacion de movimiento.
Luego se determinara la posicion del perihelio y la velocidad angular de precesion. Finalmente se ploteara el grafico de la energia en funcion
del tiempo y de la trayectoria.

'''

vy0=0.3
condicion_inicial = sp.array([10, 0, 0, vy0])
p = Planeta(condicion_inicial, 10**(-2.375))
t_final =  180*30                      #factor 30 por las orbitas
pasos = 40000
dt= t_final / (float)(pasos)

x = sp.zeros(pasos)
y = sp.zeros(pasos)
vx = sp.zeros(pasos)
vy = sp.zeros(pasos)
r = sp.zeros(pasos)
energia = sp.zeros(pasos)
perihelio = [[], [],[] ]             # se trabajara con dt,dy y dx

''' Condiciones iniciales'''
[x[0],y[0],vx[0],vy[0]] = condicion_inicial
r[0] = sp.sqrt(x[0]**2+y[0]**2)

Este script utiliza la clase Planeta para resolver la órbita de un planeta de
masa 1 kg, orbitando un cuerpo central de masa M = 1.498 x 10^10 kg. El
movimiento es determinado usando el método de Verlet. En este script se añade
el parámetro correspondiente a la corrección relativista para órbitas cercanas
al cuerpo central. Al final se presentan los resultados en forma de gráficos de
la órbita y de la energía en el tiempo.
'''

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt

#Solucion a trayectoria
condicion_inicial = [10.0, 0.0, 0.0, 0.3]

Mercury = Planeta(condicion_inicial,10**(-2.350))
x_pos = [Mercury.y_actual[0]]
y_pos = [Mercury.y_actual[1]]
energia = [Mercury.energia_total()]
dt=0.01
Mercury.avanza_rk4(dt)
for i in range(600000):
    x_pos.append(Mercury.y_actual[0])
    y_pos.append(Mercury.y_actual[1])
    energia.append(Mercury.energia_total())
    Mercury.avanza_verlet(dt)

x_pos = np.array(x_pos)
y_pos = np.array(y_pos)
energia = np.array(energia)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from planeta import Planeta
import numpy as np
import matplotlib.pyplot as plt


condicion_inicial = [10, 0, 0, 0.4]

p = Planeta(condicion_inicial)

dt = 0.1
iteraciones = 60000
x = np.zeros(iteraciones)
y = np.zeros(iteraciones)
vx = np.zeros(iteraciones)
vy = np.zeros(iteraciones)
energia = np.zeros(iteraciones)
tiempo = np.zeros(iteraciones)

x[0] = p.y_actual[0]
y[0] = p.y_actual[1]
energia[0] = p.energia_total()
tiempo[0] = p.t_actual
vx[0] = p.y_actual[2]
vy[0] = p.y_actual[3]

#Verlet necesita una iteracion extra
p.avanza_rk4(dt)
x[1] = p.y_actual[0]

# -*- coding: utf-8 -*-
'''
Script que utiliza Runge_kutta 4 para calcular y graficar aproximadamente
5 orbitas. Además grafica energia vs tiempo.
'''

import numpy as np
import matplotlib.pyplot as plt
from planeta import Planeta

condicion_inicial = [10, 0, 0, 0.3] #x, y, vx, vy

p = Planeta(condicion_inicial) #se establece condicio inicial como actuales

t_total=900 #aproximadamente 5 orbitas
# se crean listas vacias de x, y, Energia
x=[]
y=[]
Energia=[]
#Se agregan condiciones iniciales
x.append(condicion_inicial[0])
y.append(condicion_inicial[1])
Energia.append(p.energia_total())

dt=1
t=1 #ya se realizo el primer calculo
while t<=t_total:
    p.avanza_rk4(dt)
    xf, yf, vxf, vyf = p.y_actual
    x.append(xf)
    y.append(yf)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from planeta import Planeta
import matplotlib.pyplot as plt
import numpy as np

alpha = 0
condicion_inicial = [10, 0, 0, 0.2]

p = Planeta(condicion_inicial,alpha)
print "Energia Total", p.energia_total()
contador = 0
x=np.zeros(100001)
y=np.zeros(100001)
t=np.linspace(0,2000,100001)
E_total = np.zeros(100001)
while  contador<=100000:
    E_total[contador] = p.energia_total()
    x[contador] = p.y_actual[0]
    y[contador] = p.y_actual[1]
    p.avanza_verlet(0.02)
    contador +=1

print np.shape(t),np.shape(x)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Orbitas con Metodo de Verlet")
plt.plot(x,y,'-')
plt.savefig("Orbitas con Metodo de Verlet")
plt.show()