fork(1) download 
  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.         fouradded=vector.add(comebine1, comebine2)
  31.         return fouradded 
  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.         term5=self.smultiply(c5*vec1[2],(0,0,1))
  48.         firstfouradded=self.addfour(term1,term2,term3,term4)
  49.         allfiveadded=self.add(firstfouradded, term5)
  50.         out=self.smultiply(h, allfiveadded)
  51.         return  out
  52.     def rungekutta1(self,vec1,vec):
  53.         k11=self.smultiply(h,vec)
  54.         k11_half=self.smultiply(.5, k11)
  55.         y12=self.add(vec, k11_half)
  56.         k12=self.smultiply(h,y12)
  57.         k12_half=self.smultiply(.5, k12)
  58.         y13=self.add(vec, k12_half)
  59.         k13=self.smultiply(h,y13)
  60.         y14=self.add(vec, k13)
  61.         k14=self.smultiply(h,y14)
  62.         k12_double=self.smultiply(2.0, k12)
  63.         k13_double=self.smultiply(2.0, k13)
  64.         k1s=self.addfour(k11, k12_double, k13_double, k14)
  65.         n_next=self.add(vec1, self.smultiply(1.0/6, k1s))
  66.         return n_next
  67.     def rungekutta2(self, vec1, vec2):
  68.         k21=self.function(vec1, vec2)
  69.         n22=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  70.         d_n22=self.add(vec2,self.smultiply(.5,k21))
  71.         k22=self.function(n22,d_n22)
  72.         n23=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  73.         d_n23=self.add(vec2,self.smultiply(.5,k22))
  74.         k23=self.function(n23,d_n23)
  75.         n24=self.add(vec1,self.smultiply(h,(1,1,1)))
  76.         d_n24=self.add(vec2,k23)
  77.         k24=self.function(n24,d_n24)
  78.         k22_double=self.smultiply(2.0,k22)
  79.         k23_double=self.smultiply(2.0,k23)
  80.         k2s=self.addfour(k21, k22_double, k23_double, k24)
  81.         d_n_next=self.add(vec2, self.smultiply(1.0/6, k2s))
  82.         return d_n_next
  83.  
  84. def plotlines(x ,y ,z):
  85.         line = ax.plot([x,0],[y,0],[z,0],color='#0000A0',marker='.')
  86.         ax.set_xlabel('X-axis')
  87.         ax.set_ylabel('Y-axis')
  88.         ax.set_zlabel('Z-axis')
  89.         #del ax.lines[0]
  90.         matplotlib.pyplot.show()
  91.  
  92.  
  93.  
  94.  
  95. def plotdots(x,y,z):
  96.     l = ax.scatter(x, y, z, c='#387C44',marker='.')
  97.     del ax.lines[0]
  98.     ax.set_xlabel('X-axis')
  99.     ax.set_ylabel('Y-axis')
  100.     ax.set_zlabel('Z-axis')
  101.     draw()
  102.     #matplotlib.pyplot.show()
  103.  
  104.  
  105. def plotvel(dx, dy, dz):
  106.         print dx,dy,dz
  107.         ax.scatter(dx, dy, dz, c='m')
  108.         ax.set_xlabel('X-axis')
  109.         ax.set_ylabel('Y-axis')
  110.         ax.set_zlabel('Z-axis')
  111.         matplotlib.pyplot.show()
  112.  
  113. def xyproj(xs,ys):
  114.     print xs,ys
  115.     #plot(xs,ys,'o')
  116.     plotpos(xs ,ys ,0)
  117.     show()
  118.  
  119. def drawstaticplot(m,n, d_n):    
  120.     for i in range(0,m):
  121.         n=vector.rungekutta1(n, d_n)
  122.         d_n=vector.rungekutta2(n, d_n)
  123.         x1 = n[0]    
  124.         y1 = n[1]
  125.         z1 = n[2]
  126.         #print x1,y1,z1
  127.         xarray.append(x1)
  128.         yarray.append(y1)
  129.         zarray.append(z1)
  130.     for j in range(0,m-20):
  131.         if j%20 == 0:
  132.             ax.plot([xarray[j],xarray[j+20]],[yarray[j],yarray[j+20]],[zarray[j],zarray[j+20]],color='#817339',marker='.')
  133.             matplotlib.pyplot.show()
  134.  
  135.  
  136. def drawdynamicplot(m,n,d_n):
  137.     for i in range(0,m):
  138.             n=vector.rungekutta1(n, d_n)
  139.             d_n=vector.rungekutta2(n, d_n)
  140.             x1 = n[0]
  141.             y1 = n[1]
  142.             z1 = n[2]
  143.             print x1,y1,z1
  144.             if i%20 == 0:
  145.                 print "-----------------------------"
  146.                 plotlines(x1,y1,z1)
  147.                 #plotdots(x1,y1,z1)
  148.  
  149. def realtimeplotline(m,n,d_n):
  150.     for i in range(0,m):
  151.             n=vector.rungekutta1(n, d_n)
  152.             d_n=vector.rungekutta2(n, d_n)
  153.             x1 = n[0]
  154.             y1 = n[1]
  155.             z1 = n[2]
  156.             print x1,y1,z1
  157.             if i%20 == 0: 
  158.                 print "--------------------mark------------------"
  159.                 index = ((i/20))%2
  160.                 #print index
  161.                 if index ==0:
  162.                     x2 = n[0]
  163.                     y2 = n[1]
  164.                     z2 = n[2]
  165.                 if (index) == 1:
  166.                     x3 = n[0]
  167.                     y3 = n[1]
  168.                     z3 = n[2]
  169.                     ax.plot([x2,x3],[y2,y3],[z2,z3],color='#817339',marker='.')
  170.                     matplotlib.pyplot.show()
  171.  
  172.  
  173.  
  174. fig = plt.figure()
  175. ax = fig.add_subplot(111, projection='3d')
  176.  
  177. theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
  178.  
  179. c1=-1.0
  180. c2=-1.0
  181. c3=1.0
  182. c4=-1.0
  183. c5=0
  184. h=.004
  185. vector=vec(1,2,3)
  186. allone=vec(1.0,1.0,1.0)
  187. n=(1.0,1.0,1.0)
  188. d_n=(1.0,-1.0,0.0)
  189. ion()
  190. x = n[0]    
  191. y = n[1]
  192. z = n[2]
  193. n=(1.0,1.0,1.0)
  194. d_n=(1.0,-1.0,0.0)
  195.  
  196.  
  197. m=10000
  198. xarray=[]
  199. yarray=[]
  200. zarray=[]
  201.  
  202. #drawstaticplot(m,n, d_n)
  203. drawdynamicplot(m,n, d_n)
  204. #realtimeplotline(m,n, d_n)
  205.  
Runtime error #stdin #stdout 0.09s 10920KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/0odwC
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!