VPython

[Miro]

ora, ja agora fica aqui o codigo da simples apresentação do nosso projecto para a fisica@uc, mas isto num windows é feio, usem um mac!!!! no windows texturas nao aparecem

o sistema solar e mesmo pa enfeitar, nao liguem a tonelada de imperfeiçoes, o objecto que roda e fixe para ver as difrenças no brilho dependendo da nossa posição que era todo o nosso objectivo

[Tharis]

Nice!

[Miro]

alguem me consegue por isto a dar? não consigo ver qual é o erro..... estava a dar....e deixou de dar.....sem razao aparente

from __future__ import division
from visual import *


def make_time_string(t):
"Accept a number of seconds, return a relative string."
if t < 60: return "%02i seconds"%t
mins,secs = divmod(t, 60)
if mins < 60: return "%02i minutes %02i seconds"%(mins,secs)
hours, mins = divmod(mins,60)
if hours < 24: return "%02i hours %02i minutes"%(hours,mins)
days,hours = divmod(hours,24)
return "%02i days %02i hours"%(days,hours)


#window
scene.center = (0,0,0)
scene.width = 800
scene.height = 600
scene.range = (1.0e9,1.0e9,1.0e9)

#constantes
dt = input('dt ')
totalseconds = 0
G = 6.67e-11
K = 8.9875517873681764e9
electrao = -1.602e-19
protao = -1.602e-19
#objectos

earth = sphere(pos=(0,0,0))

terra =6.38e6
lua =1.74e6
sol = 6.955e8

earth.radius=input('raio corpo1 ')
earth.color=color.blue

terra =5.98e24
lua =7.36e22
sol =1.9891e30

earth.mass = input ('massa corpo1 ')
earth.carga = input ('carga corpo1 (electrao / protao) ')

terralua = 3.84e8
terrasol = 1.496e11

afastamento = input ('distancia entre 2 corpos ')

lua = sphere(pos=(0,afastamento,0))

terra =6.38e6
lua =1.74e6
sol = 6.955e8

lua.radius = input ('raio corpo2 ')
lua.color = color.white

terra =5.98e24
lua =7.36e22
sol =1.9891e30

lua.mass = input ('massa corpo 2 ')
lua.carga = input ('carga corpo2 (electrao / protao) ')

info = label(pos=(0,afastamento + 100000000,0))


#velocidade inicial da lua
initialAcceleration = ((G * earth.mass * lua.mass / 3.84e8**2)+(K*earth.carga*lua.carga/3.84e8**2)) / lua.mass

calculatedVelocity = (initialAcceleration * 3.84e8)**.5
velocidadedalua = calculatedVelocity

lua.velocity = velocidadedalua * vector(1,0,0)
earth.velocity = vector (0,0,0)
#main

finished = False
while not finished:

totalseconds += dt
rate(10000)

distancia = mag(earth.pos - lua.pos)

ForceGravityDir = norm(earth.pos - lua.pos)

ForceGravityMag = G * (earth.mass * lua.mass) / distancia**2
ForceElectricMag = K * (earth.carga * lua.carga) / distancia**2


ForceGravity = ForceGravityMag * ForceGravityDir
ForceGravityT = -1 * ForceGravityMag * ForceGravityDir #contrario a do outro sat leva -1


if lua.carga * earth.carga < 1:
ForceElectric = ForceElectricMag * ForceGravityDir *-1
ForceElectricT = ForceElectricMag * ForceGravityDir
if lua.carga * earth.carga > 1:
ForceElectric = ForceElectricMag * ForceGravityDir
ForceElectricT = ForceElectricMag * ForceGravityDir *-1

acelE = ForceElectric /lua.mass
acelG = ForceGravity /lua.mass


lua.acceleration = (ForceGravity+ForceElectric) / lua.mass
lua.velocity += lua.acceleration * dt
lua.pos += lua.velocity * dt + .5 * lua.acceleration * dt**2

earth.acceleration = (ForceGravityT+ForceElectricT) / earth.mass
earth.velocity += earth.acceleration * dt
earth.pos += earth.velocity * dt + .5 * earth.acceleration * dt**2


message = "tretas para dizer ao user, tenho de por isto comprido para a janela nao bater mal"
message += "\ndistancia " + str(((distancia-earth.radius)/1000))
message += "\nvelocidade " + str (mag(lua.velocity))
message += "\naceleracao total " + str (mag(lua.acceleration))
message += "\naceleracao gravi " + str (acelG)
message += "\naceleracao electri " + str (acelE)
message += "\nforca total " + str (mag(ForceGravity+ForceElectric))
message += "\nforca electrica " + str (mag(ForceElectric))
message += "\nforca gravitica " + str (mag(ForceGravity))
message += "\nTempo: " + make_time_string(totalseconds)
info.text = message

[jap]

O ficheiro não tem formatação (identação)?

[Bruno Oliveira]

Estou agora a mexer no Vpython, e gostava de criar uma animação simples, de um edíficio de 1 piso, com 4 pilares a oscilar sem estar sujeito à acção de forças externas.

Estou a usar 5 objectos, box, para simularem os pilares e a laje de piso.

Já consegui simular a oscilação da laje de piso resolvendo a eq. diferencial que rege o seu movimento oscilatório, e gostava de saber como posso associar à laje, alguns dos pontos das caixas que simulam os pilares... Será possível?

[ajoanagil]

Esta linguagem é a mesma utilizada nos arduinos?

[Tharis]

Não. A linguagem usada nos Arduinos é uma "espécie" de C simplificado.

[ruifm]

Fiz esta coisinha como uma preparação para o hyper-mega projecto que hei-de fazer conjuntamente com people daqui.

O código tá um bocadinho forçado e tá assim à papo-seco as coisas lá.... Utilizei a 1a lei de kepler e pouco mais de interessante...

Código: Seleccionar Todos

#!/usr/bin/env python2.6.2
#
#  Solar.py v0.1 (following the rule "if it doesn't work, name it v0.1")
# 
#
#  Created by Francisco Huhn on 6/8/09
#  Should mention Satelite.py by Joao Casimiro
#  Copyright (c) 2009. All rights reserved.
#

""" This is a script that simulates the solar system planet's
   orbits. It takes the first of Kepler's laws of planetary
   motion which calculates the orbit's radius in a given moment
   for a given angle. If you don't know much of the subject
   I suggest you read and calculate wikipedia's page on the
   subject: http://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion """


# TODO:
#
# - An abstract way of putting bodies in scene
# - Clean that mess down there (aka find a pythonic and clean way to do the job)
# - Planet sattelite support
# - A decent stat shower (in rt)
# - Clean the code up
# - User interaction
# - Maybe a shell and stuff
# - And lots of other stuff


######## MODULES ########

import visual as v
from math import *
from time import sleep
from thread import start_new_thread as nthread

#########################

####### CONSTANTS #######

G = 6.67e-11 # Univertial Gravitational Constant
pi = 3.141592653589793238462643 # pi | I could import it from math, but ...

#########################

####### SCENE CFG #######

v.scene.title = "Solar System"
v.scene.width = 1200
v.scene.height = 600

#########################

#def shell():
#   global planets, camera
#   
#   while True:
#      cmd = raw_input()
#      try:
#         cmd, target = cmd.split()
#      except:
#         pass
#   
#      if cmd == "center":
#         for planet in planets:
#            if planet == target:
#               camera = planet
#         if target == "sun":
#            camera = target
      

if __name__ == "__main__":
   
   planets = ["mercury", "venus", "earth", "mars", "jupiter"]
   
   
   """ data is a dictionary of lists of the following form [(semi-major axis, 0, 0), radius, color, aphelion, eccentricity, period (in days)]
      0 = pos
      1 = radius
      2 = color
      3 = aphelion
      4 = eccentricity
      5 = period
   """
   
   data = {"mercury":
      {"sm-axis": 57909100, "radius": 2439.7, "color": v.color.red, "aphelion": 69816900, \
      "eccentricity": 0.205630, "tperiod": 87.969, "rperiod": 58.646},
      
         "venus":
      {"sm-axis": 108208930, "radius": 6501.8, "color": v.color.white, "aphelion": 108942109, \
      "eccentricity": 0.0068, "tperiod": 224.70069, "rperiod": 243.0185},
      
         "earth":
      {"sm-axis": 147098074, "radius": 6371, "color": v.color.blue, "aphelion": 149597887.5, \
      "eccentricity": 0.016710219, "tperiod": 365.256366, "rperiod": 0.99726968 },
      
         "mars":
      {"sm-axis": 227939100, "radius": 3396, "color": v.color.orange, "aphelion": 249209300, \
      "eccentricity": 0.093315, "tperiod": 686.971, "rperiod": 1.025957},
      
         "jupiter":
      {"sm-axis": 778547200, "radius": 71492 / 10., "color": v.color.yellow, "aphelion": 816520800, \
      "eccentricity": 0.048775, "tperiod": 4311.572, "rperiod": 9.925 / 24} }

#         "moon":
#      { "sm-axis": 384399, "radius": 1737, "color": (0.5,0.5,0.5) , "aphelion": 405696, "eccentricity": 0.054, "tperiod": 27.321582 } }


   """ ps is a list with the values of p which is the semi-lactus rectum (what a weird name) of the ellipses """
   
   ps = {}
   for planet in data.keys():
      ps[planet] = data[planet]["aphelion"] * (1 - data[planet]["eccentricity"]**2)

      
   """ bodies is a list with vpython's spheres (the planets) """
   
   bodies = {}   
   for planet in data.keys():
      bodies[planet] = v.sphere(pos = (data[planet]["aphelion"],0,0) , radius = data[planet]["radius"] * 500, color = data[planet]["color"])


   """ trails is a list with vpython's curves (trails of the planets) """

   trails = {}
   for planet in data.keys():
#      trails[planet] = v.curve(color = data[planet]["color"])
      trails[planet] = v.curve(color = (1,1,1))
      
   """ tetas is a list with the values of teta which is the angle between the aphelion and r, the radius vector """
      
   tetas = {}
   alphas = {}
   for planet in data.keys():
      alphas[planet] = 0
      tetas[planet] = 0
   
   
   """ sun sphere object from vpython """
   sun = v.sphere(pos = (0,0,0), radius = 109 * 6371 * 10, color = v.color.red)
   
   
   """ xx and yy axis """
   xx = v.curve( x = (-800000000, 800000000) )
   yy = v.curve( y = (-800000000, 800000000) )
   v.label(pos = (800000000,0,0), text = "x")
   v.label(pos = (0,800000000,0), text = "y")
      
##   """ a new frame with the data of each planet """
##   data_frame = v.display(title = "Stats", width = 600, height = 200)
##   data_text = ""   
##
##   ############## LAME NOT SO PYTHONIC CODE ##############
##   # When I find a pythonic way to do it, I'll change it #
##   #######################################################
##
##   words = ["planet", "sm-axis", "radius", "aphelion", "eccentricity", "period (in days)"]
##
##   datas = [planets, [data[k][0][0] for k in data.keys()]]
##   for i in [1,3,4,5]:
##      datas.append( [ data[k][i] for k in data.keys() ] )
##   
##   for i in range(len(words)):
##      data_text += words[i] + "     "
##   data_text += "\n"
##      
##   for i in range(len(datas[0])):
##      for j in range(len(datas)):
##         data_text += str(datas[j][i]) + "      "
##      data_text += "\n"
##   
##   #######################################################
##   #######################################################
##   
##   v.label( display = data_frame, text = data_text )

   # Funny :p
   bodies["earth"].color = (1,1,1)
   bodies["earth"].material = v.materials.earth
   bodies["earth"].rotate(angle = pi/2, axis = (1,0,0))
   
   while True:
      v.rate(25)
      v.scene.center = bodies["earth"].pos
      
      for k in data.keys():
      
         bodies[k].rotate(angle = 2*pi / data[k]["rperiod"], axis = (0,0,1))#, axis = (0,1,0), origin = bodies[k].pos)
         tetas[k] += (2*pi) / data[k]["tperiod"] # add angular displacement
         
         r = ps[k] / (1 + data[k]["eccentricity"] * cos(tetas[k])) # radius
         

         bodies[k].x = cos(tetas[k]) * r
         bodies[k].y = sin(tetas[k]) * r
         
#         if tetas[k] >= 2*pi: # clean last orbit
#            tetas[k] = tetas[k] % (2*pi)
#            trails[k].pos = [(data[k]["aphelion"],0,0)]
         
         trails[k].append(pos = bodies[k].pos) # orbit's trail

eu tambem tive a ideia de fazer um sistema solar em vpython, e ja o consegui usando apenas a 1ª Lei de Kepler e dados astronomicos que se encontram na net.

Mas achei mais interessante criar um programa que usasse APENAS as Leis de Newton para criar o sistema solar, e assim confirmar que as orbitas sao elipticas como Kepler afirmou. E consegui, em pequena escala, com 2 esferas. Se desse uma determinada velocidade o corpo atingia um MCU, caso contrario a orbita seria uma elipse com o seu eixo constantemente a variar. Se a velocidade fosse muito baixa, como é de esperar o corpo iria de encontro ao corpo de maior massa.
Como tudo funcionou na perfeiçao na micro escala, resolvi colocar valores reais. Eu tenho consciencia que desta forma seria dificil ver a Terra em comparaçao com o Sol, e que o movimento de Translacao seria praticamente invisivel porque demora 1 ano xD. Para adicionar um pouco mais de realismo, criei ainda movimento de rotaçao na Terra, com esta declinada 23.4º.
O movimento começa no afelio. Por essa razao achei interessante obter dados relativamente a diminuicao da distancia, e ao aumento de velocidade. Por isso vao aparecer 2 graficos espetaculares em pop up.
Espero que gostem!!

Código: Seleccionar Todos

from visual.graph import * #importar graficos tambem
import math
#este programa permite exprimir o movimento da Terra dado o seu raio, a sua distancia ao sol, e a as massas dos mesmos. a partir dai e possivel controlar e prever o movimento de qualquer outro corpor com qualquer outra velocidade

#janela
scene.autoscale=0 #tirar o zoom irritante
scene.title='Earth orbit by Newton'
scene.width=1680 #resolucao do meu PC
scene.height=1050
scene.range=(2*10**11,2*10**11,2*10**11)
scene.lights=[] #importante...desligar luzes iniciais
scene.show_rendertime=True #mostrar velocidade dos ciclos (opcional, para controlo)

#corpos
scene.select()
sun=sphere(pos=(0,0,0), radius=6.96*10**8, color=color.yellow, mass=1.99*10**30, material=materials.emissive)
earth=ellipsoid(pos=(0,0,152098233000), size=(2*6378137,2*6356752.3,2*6378137), mass=5.98*10**24, material=materials.earth, make_trail=True, trail_type="points", interval=100, retain=999999)
lamp=local_light(pos=sun.pos, color=color.white) #luz vinda do sol
earth.trail_object.color=color.blue #deixar rasto da terra azul em vez de branco
#graficos
gd1=gdisplay(title='Distance(t)', xtitle='time(s)', ytitle='distance(m)') # 1a janela de grafico
f1=gcurve(color=color.red, display=gd1.display)
value1=label(display=gd1.display) #mostrador de valor
gd2=gdisplay(title='Speed(t)', xtitle='time(s)', ytitle='speed(m/s)') # 2a janela de grafico
f2=gcurve(color=color.green, display=gd2.display)
value2=label(display=gd2.display) #mostrador de valor

#constantes
G=6.67*10**(-11)
sun.vel=vector(0,0,0)
earth.vel=vector(29291,0,) # velocidade inicial no afelio
earth.rotate(angle=0.408407045, axis=(0,0,1)) # declinacao da Terra
t=0.
speed='Earth speed: ' #para mostradores de resultados... soon
distance='Distance: '

#motion
while True:
    rate(100)
    scene.center=earth.pos #melhor perspectiva
    d=math.sqrt((earth.pos.x-sun.pos.x)**2 + (earth.pos.y-sun.pos.y)**2 + (earth.pos.z-sun.pos.z)**2) #distancia sol-terra
    F=(G*sun.mass*earth.mass)/(d**2) #lei da gravitacao universal
    a=F/earth.mass # 2a lei de newton
    a2=F/sun.mass # para o sol
    accelaration=vector(sun.pos.x-earth.pos.x,sun.pos.y-earth.pos.y,sun.pos.z-earth.pos.z) #vector direcional da normal
    k=a/math.sqrt((accelaration.x**2)+(accelaration.y**2)+(accelaration.z**2)) # k que ajusta o modulo da aceleracao
    centripeta=vector(accelaration.x*k,accelaration.y*k,accelaration.z*k) # ajuste do modulo para ter a verdadeira aceleracao centripeta
    earth.vel+=centripeta
    earth.pos+=earth.vel
    earth.rotate(angle=((math.pi*2)/(24*360000)), axis=(0,1,0)) # movimento de rotacao da Terra, so visivel se o rate for aumentado, ou mesmo tirado
    accelaration2=vector(earth.pos.x-sun.pos.x,earth.pos.y-sun.pos.y,earth.pos.z-sun.pos.z) # repetir tudo para o sol (nao e desprezavel)
    k2=a2/math.sqrt((accelaration2.x**2)+(accelaration2.y**2)+(accelaration2.z**2))
    centripeta2=vector(accelaration2.x*k2,accelaration2.y*k2,accelaration2.z*k2)
    sun.vel+=centripeta2
    sun.pos+=sun.vel
    lamp.pos=sun.pos
    f1.plot(pos=(t,d-152098233000)) #grafico distancia-tempo
    value1.text=distance+str(d)
    value1.pos=(t,d-152098233000)
    f2.plot(pos=(t,mag(earth.vel)-29291)) #grafico velocidade-tempo
    value2.text=speed+str(mag(earth.vel))
    value2.pos=(t,mag(earth.vel)-29291)
    t+=0.01 #gerador

só há um problema. O modo visual fica muito lento passado cerca de 300 segundos e eu não uso recursao nem nada. Precisava de uma granda ajuda nesta parte, pelos experts da memoria

[filipematos]

Mas qual foi a ideia de a terra ser maior do que o SOL?

[ruifm]

Mas qual foi a ideia de a terra ser maior do que o SOL?

como assim?

[filipematos]

Estou a analisar o teu trabalho geocentrista

[filipematos]

Porque deste uma forma estranha à terra? É melhor pores uma esfera não?

[ruifm]

Estou a analisar o teu trabalho geocentrista

não está geocentrista!
é puro newton...
o ponto de vista é que é da Terra, para ser mais familiar.

A Terra não está com uma forma estranha, está 'oval' (o tal ellipsoid) que é a forma mais realista. Tem um raio equatorial diferente do raio polar.
eu quero é mesmo perceber porque razao o modo visual fica lento após uns minutos

[Bruno Oliveira]

Experimenta alterar o teu argumento do rate, para ver que tipo de alteração ele provoca na simulação

Ou então pode estar relacionado com o t. Tu começas em 0, e incrementas 0.01 por iteração, mas quando é que acaba a contagem do tempo? Por outras palavras, não faria mais sentido correres a simulação durante um tempo t previamente fixado? Assim o teu while true mudaria para while t < tmax...

Isto foi o que me saltou mais à vista ao ler o teu código na diagonal... Espero ter tempo de o ler a sério, quando o puder fazer, avisarei!

Bruno

[ruifm]

Experimenta alterar o teu argumento do rate, para ver que tipo de alteração ele provoca na simulação

Ou então pode estar relacionado com o t. Tu começas em 0, e incrementas 0.01 por iteração, mas quando é que acaba a contagem do tempo? Por outras palavras, não faria mais sentido correres a simulação durante um tempo t previamente fixado? Assim o teu while true mudaria para while t < tmax...

Isto foi o que me saltou mais à vista ao ler o teu código na diagonal... Espero ter tempo de o ler a sério, quando o puder fazer, avisarei!

Bruno

o programa corre bem nos primeiros minutos. Se limitasse o tempo eu tenho a certeza que o problema desaparecia...
Mas após uns 10 minutos o render time, que tenho a ideia de ser o intrevalo de tempo entre cada ciclo, é muito grande (500 ms) comparativamente com o mesmo no inicio (30 ms, que é o minimo).
mas eu queria mesmo que o programa durasse o tempo que me apetecesse. A variavel t so existe mesmo para o grafico, eu poderia retira la e corria tudo bem.
Terera alguma coisa a ver com o trail? Se eu o retirar ou colocar um intrevalo maior...
Sera que é da textura da Terra que ocupa memoria?
Ou é a quantidade de movimentos e interaçoes que são demais para um ciclo? E se assim for, porque é que tudo corre bem no inicio?