fork(5) download
copy 
  1. import math
  2. import matplotlib.pyplot
  3. import matplotlib.pyplot as plt
  4. import matplotlib as mpl
  5. from mpl_toolkits.mplot3d import Axes3D
  6. from matplotlib.collections import PolyCollection
  7. from matplotlib.colors import colorConverter
  8. import numpy as np
  9. import matplotlib.pyplot as plt
  10. from numpy import *
  11. from pylab import *
  12. import pylab
  13.  
  14.  
  15. class vec(object):
  16.     def __init__(self, x, y, z):
  17.         self.vec=(x, y, z)        
  18.  
  19.     def unit(self, vec):
  20.         dist=math.sqrt(vec[0]**2+vec[1]**2+vec[2]**2)
  21.         u=(vec[0]/dist,vec[1]/dist,vec[2]/dist)
  22.         return u
  23.  
  24.     def add(self, vector1, vector2):
  25.         sum=(vector1[0]+vector2[0], vector1[1]+vector2[1], vector1[2]+vector2[2])
  26.         return sum
  27.     def addfour(self, vec1, vec2, vec3, vec4):
  28.         comebine1=vector.add(vec1,vec2)
  29.         comebine2=vector.add(vec3,vec4)
  30.         outcome=vector.add(comebine1, comebine2)
  31.         return outcome   
  32.     def smultiply(self,scalar, vec):
  33.         product=(scalar*vec[0],scalar*vec[1],scalar*vec[2])
  34.         return product
  35.     def innerproduct(self, vec1, vec2):
  36.         result=(vec1[0]*vec2[0]+vec1[1]*vec2[1]+vec1[2]*vec2[2])
  37.         return result
  38.     def crossproduct(self, vec1, vec2):
  39.         res=(vec1[1]*vec2[2]-vec1[2]*vec2[1], vec1[2]*vec2[0]-vec1[0]*vec2[2], vec1[0]*vec2[1]-vec1[1]*vec2[0])
  40.         return res
  41.     def function(self, vec1, vec2):
  42.         coef1=c1*inner(vec2, vec2)
  43.         term1=self.smultiply(coef1, vec2)
  44.         term2=self.smultiply(c2,self.crossproduct(vec2,vec1))
  45.         term3=self.smultiply(c3,(0,0,1))
  46.         term4=self.smultiply(c4*vec1[2],vec1)
  47.         sumoffour=self.addfour(term1,term2,term3,term4)
  48.         out=self.smultiply(h, sumoffour)
  49.         return  out
  50.     def rungekutta1(self,vec1,vec):
  51.         k11=self.smultiply(h,vec)
  52.         k11_half=self.smultiply(.5, k11)
  53.         y12=self.add(vec, k11_half)
  54.         k12=self.smultiply(h,y12)
  55.         k12_half=self.smultiply(.5, k12)
  56.         y13=self.add(vec, k12_half)
  57.         k13=self.smultiply(h,y13)
  58.         y14=self.add(vec, k13)
  59.         k14=self.smultiply(h,y14)
  60.         k12_double=self.smultiply(2.0, k12)
  61.         k13_double=self.smultiply(2.0, k13)
  62.         k1s=self.addfour(k11, k12_double, k13_double, k14)
  63.         n_next=self.add(vec1, self.smultiply(1.0/6, k1s))
  64.         return n_next
  65.     def rungekutta2(self, vec1, vec2):
  66.         k21=self.function(vec1, vec2)
  67.         n22=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  68.         d_n22=self.add(vec2,self.smultiply(.5,k21))
  69.         k22=self.function(n22,d_n22)
  70.         n23=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  71.         d_n23=self.add(vec2,self.smultiply(.5,k22))
  72.         k23=self.function(n23,d_n23)
  73.         n24=self.add(vec1,self.smultiply(h,(1,1,1)))
  74.         d_n24=self.add(vec2,k23)
  75.         k24=self.function(n24,d_n24)
  76.         k22_double=self.smultiply(2.0,k22)
  77.         k23_double=self.smultiply(2.0,k23)
  78.         k2s=self.addfour(k21, k22_double, k23_double, k24)
  79.         d_n_next=self.add(vec2, self.smultiply(1.0/6, k2s))
  80.         return d_n_next
  81.  
  82. def plotlines(x ,y ,z):
  83.         line = ax.plot([x,0],[y,0],[z,0],color='#0000A0',marker='.')
  84.         #ax.lines.pop(0)
  85.         ax.set_xlabel('X-axis')
  86.         ax.set_ylabel('Y-axis')
  87.         ax.set_zlabel('Z-axis')
  88.         #del ax.lines[0]
  89.         matplotlib.pyplot.show()
  90.  
  91.  
  92.  
  93.  
  94. def plotdots(x,y,z):
  95.     l = ax.scatter(x, y, z, c='#387C44',marker='.')
  96.     ax.set_xlabel('X-axis')
  97.     ax.set_ylabel('Y-axis')
  98.     ax.set_zlabel('Z-axis')
  99.     matplotlib.pyplot.show()
  100.  
  101.  
  102. def plotvel(dx, dy, dz):
  103.         print dx,dy,dz
  104.         ax.scatter(dx, dy, dz, c='m')
  105.         ax.set_xlabel('X-axis')
  106.         ax.set_ylabel('Y-axis')
  107.         ax.set_zlabel('Z-axis')
  108.         matplotlib.pyplot.show()
  109.  
  110. def xyproj(xs,ys):
  111.     print xs,ys
  112.     #plot(xs,ys,'o')
  113.     plotpos(xs ,ys ,0)
  114.     show()
  115.  
  116. def drawstaticplot(m,n, d_n):    
  117.     for i in range(0,m):
  118.         n=vector.rungekutta1(n, d_n)
  119.         d_n=vector.rungekutta2(n, d_n)
  120.         x1 = n[0]    
  121.         y1 = n[1]
  122.         z1 = n[2]
  123.         #print x1,y1,z1
  124.         xarray.append(x1)
  125.         yarray.append(y1)
  126.         zarray.append(z1)
  127.     for j in range(0,m-20):
  128.         if j%20 == 0:
  129.             ax.plot([xarray[j],xarray[j+20]],[yarray[j],yarray[j+20]],[zarray[j],zarray[j+20]],color='#817339',marker='.')
  130.             matplotlib.pyplot.show()
  131. def connectpoints(x1,y1,z1,x2,y2,z2):
  132.     ax.plot([x2,x2],[y2,y2],[z2,z2],color='#817339',marker='.')
  133.     matplotlib.pyplot.show()
  134.  
  135.  
  136. def realtimeplotline(m,n,d_n):
  137.     x1=n[0]
  138.     y1=n[1]
  139.     z1=n[2]
  140.     x2=0
  141.     y2=0
  142.     z2=0
  143.     for i in range(0,m):
  144.             n=vector.rungekutta1(n, d_n)
  145.             d_n=vector.rungekutta2(n, d_n)
  146.             x2= n[0]
  147.             y2= n[1]
  148.             z2= n[2]
  149.             if i%20 == 0:
  150.                 connectpoints(x2,y2,z2,x1,y1,z1)
  151.                 x1=x2
  152.                 y1=y2
  153.                 z1=z2
  154.  
  155.  
  156.  
  157.  
  158.  
  159.  
  160.  
  161.  
  162.  
  163. fig = plt.figure()
  164. ax = fig.add_subplot(111, projection='3d')
  165.  
  166. theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
  167.  
  168. c1=-1.0
  169. c2=-1.0
  170. c3=1.0
  171. c4=-1.0
  172. h=.004
  173. vector=vec(1,2,3)
  174. allone=vec(1.0,1.0,1.0)
  175. n=(1.0,1.0,1.0)
  176. d_n=(1.0,-1.0,0.0)
  177. ion()
  178. x = n[0]    
  179. y = n[1]
  180. z = n[2]
  181. n=(1.0,1.0,1.0)
  182. d_n=(1.0,-1.0,0.0)
  183.  
  184.  
  185. m=10000
  186. xarray=[]
  187. yarray=[]
  188. zarray=[]
  189.  
  190. #drawstaticplot(m,n, d_n)
  191. realtimeplotline(m,n, d_n)
  192.  
  193.  
Runtime error #stdin #stdout 0.09s 10840KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/lQgmS
language:
Python (cpython 2.7.16)
created:
13 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!