fork download
copy 
  1. import math
  2. #import matplotlib.pyplot
  3. #import matplotlib as mpl
  4. from mpl_toolkits.mplot3d import Axes3D
  5. from matplotlib.collections import PolyCollection
  6. from matplotlib.colors import colorConverter
  7. from matplotlib.patches import FancyArrowPatch
  8. from mpl_toolkits.mplot3d import proj3d
  9. import mpl_toolkits.mplot3d.art3d as art3d
  10. #import numpy as np
  11. import matplotlib.pyplot as plt
  12. #from numpy import *
  13. from pylab import *
  14. #import pylab
  15.  
  16.  
  17. class vec(object):
  18.     def __init__(self, x, y, z):
  19.         self.vec=(x, y, z)
  20.  
  21.     def unit(self, vec):
  22.         dist=math.sqrt(vec[0]**2+vec[1]**2+vec[2]**2)
  23.         unitvec=(vec[0]/dist,vec[1]/dist,vec[2]/dist)
  24.         return unitvec
  25.  
  26.     def add(self, vector1, vector2):
  27.         sum=(vector1[0]+vector2[0], vector1[1]+vector2[1], vector1[2]+vector2[2])
  28.         return sum
  29.  
  30.     def addfour(self, vec1, vec2, vec3, vec4):
  31.         comebine1=self.add(vec1,vec2)
  32.         comebine2=self.add(vec3,vec4)
  33.         outcome=self.add(comebine1, comebine2)
  34.         return outcome
  35.  
  36.     def smultiply(self,scalar, vec):
  37.         product=(scalar*vec[0],scalar*vec[1],scalar*vec[2])
  38.         return product
  39.  
  40.     def innerproduct(self, vec1, vec2):
  41.         result=(vec1[0]*vec2[0]+vec1[1]*vec2[1]+vec1[2]*vec2[2])
  42.         return result
  43.  
  44.     def crossproduct(self, vec1, vec2):
  45.         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])
  46.         return res
  47.  
  48.     def function(self, vec1, vec2, i):
  49.         coef1=c1*inner(vec2, vec2)
  50.         term1=self.smultiply(coef1, vec1)
  51.         term2=self.smultiply(c2,self.crossproduct(vec2,vec1))
  52.         term3=self.smultiply(c3*(1+c5*math.sin(5000*i*h)),(0,0,1))
  53.         term4=self.smultiply(c4*(1+c5*math.sin(5000*i*h))*vec1[2],vec1)
  54.         sumoffour=self.addfour(term1,term2,term3,term4)
  55.         out=self.smultiply(h, sumoffour)
  56.         return out
  57.     def rungekutta1(self,vec1,vec):
  58.         k11=self.smultiply(h,vec)
  59.         k11_half=self.smultiply(.5, k11)
  60.         y12=self.add(vec, k11_half)
  61.         k12=self.smultiply(h,y12)
  62.         k12_half=self.smultiply(.5, k12)
  63.         y13=self.add(vec, k12_half)
  64.         k13=self.smultiply(h,y13)
  65.         y14=self.add(vec, k13)
  66.         k14=self.smultiply(h,y14)
  67.         k12_double=self.smultiply(2.0, k12)
  68.         k13_double=self.smultiply(2.0, k13)
  69.         k1s=self.addfour(k11, k12_double, k13_double, k14)
  70.         n_ununit=self.add(vec1, self.smultiply(1.0/6, k1s))
  71.         n_next=self.unit(n_ununit)
  72.         return n_next
  73.     def rungekutta2(self, vec1, vec2, i):
  74.         k21=self.function(vec1, vec2, i)
  75.         n22=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  76.         d_n22=self.add(vec2,self.smultiply(.5,k21))
  77.         k22=self.function(n22,d_n22, i)
  78.         n23=self.add(vec1,self.smultiply(.5*h,(1,1,1)))
  79.         d_n23=self.add(vec2,self.smultiply(.5,k22))
  80.         k23=self.function(n23,d_n23, i)
  81.         n24=self.add(vec1,self.smultiply(h,(1,1,1)))
  82.         d_n24=self.add(vec2,k23)
  83.         k24=self.function(n24,d_n24, i)
  84.         k22_double=self.smultiply(2.0,k22)
  85.         k23_double=self.smultiply(2.0,k23)
  86.         k2s=self.addfour(k21, k22_double, k23_double, k24)
  87.         d_n_ununit=self.add(vec2, self.smultiply(1.0/6, k2s))
  88.         d_n_next=self.unit(d_n_ununit)
  89.         return d_n_next
  90.  
  91. class Arrow3D(FancyArrowPatch):
  92.     def __init__(self, xs, ys, zs, *args, **kwargs):
  93.         FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
  94.         self._verts3d = xs, ys, zs
  95.  
  96.     def draw(self, renderer):
  97.         xs3d, ys3d, zs3d = self._verts3d
  98.         xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
  99.         self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
  100.         FancyArrowPatch.draw(self, renderer)
  101.  
  102. def drawbody(x, y, z):
  103.         u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
  104.         x1=.4*np.cos(u)*np.sin(v)
  105.         y1=.4*np.sin(u)*np.sin(v)
  106.         z1=.1*np.cos(v)
  107.         theta = math.acos(math.sqrt(y**2+ z**2)/math.sqrt(x**2 + y**2 +z**2))
  108.         x2 = x1*math.cos(theta) + z1*math.sin(theta)
  109.         y2 = y1
  110.         z2 = -x1*math.sin(theta) + z1*math.cos(theta)
  111.         if theta <0:
  112.             theta1 = math.acos(z/math.sqrt(y**2+z**2))
  113.         else:
  114.             theta1=-math.acos(z/math.sqrt(y**2+z**2))
  115.             x3 = x2
  116.             y3 = y2*math.cos(theta1) - z2*math.sin(theta1)
  117.             z3 = y2*math.sin(theta1) + z2*math.cos(theta1)
  118.         ax.plot_wireframe(x3+x/2, y3+y/2, z3+z/2, color="#254117")           
  119.         ax.scatter([x3+x/2],[y3+y/2],[0], marker='_', color='black')
  120.         matplotlib.pyplot.show()
  121.  
  122. def drawgroundandlabelaxis(z, n_o):
  123.         rectangle = Rectangle((-1,-1),2,2, color='grey')
  124.         ax.add_patch(rectangle)
  125.         art3d.pathpatch_2d_to_3d(rectangle, z=z- n_o)
  126.         ax.set_xlabel('X-axis')
  127.         ax.set_ylabel('Y-axis')
  128.         ax.set_zlabel('Z-axis')     
  129.         matplotlib.pyplot.show()
  130.  
  131. def plotvel(dx, dy, dz, x, y, z):
  132.         a = Arrow3D([x,dx+x],[y,dy+y],[z,dz+z], mutation_scale=20, lw=1, arrowstyle="-|>", color='r')
  133.         ax.add_artist(a)
  134.         plt.show()
  135.         #ax.scatter(dx, dy, dz, c='#330099')
  136.         matplotlib.pyplot.show()
  137.  
  138. def plotpos(x ,y ,z, n_o, w):
  139.         ax.cla()
  140.         #drawbody(x,y,z)
  141.         line = ax.plot([x,0],[y,0],[z,z- n_o],color='#000066', marker= 'o')
  142.         line = ax.plot([0,0],[0,0],[1.1,-.5],color='g')
  143.         ax.scatter([0,0],[0,0],[0,1], marker='8', color='g')
  144.         drawgroundandlabelaxis(z, n_o)
  145.         matplotlib.pyplot.show()
  146.  
  147. def plotthetaphi(x,y,z):
  148.     phi= math.acos(z/math.sqrt(x**2+y**2+z**2))
  149.     theta = math.acos(x/math.sqrt(x**2 + y**2))
  150.     plot(theta, phi,'.',color='red')
  151.  
  152.  
  153.  
  154. def drawstaticplot(m, n, d_n, n_o):
  155.     counter=0
  156.     for i in range(0,m):
  157.         n=vector.rungekutta1(n, d_n)
  158.         d_n=vector.rungekutta2(n, d_n, i)
  159.         x1 = n[0]
  160.         y1 = n[1]
  161.         z1 = n[2]
  162.         if i%20==0:
  163.             xarray.append(x1)
  164.             yarray.append(y1)
  165.             zarray.append(z1)
  166.         x_m = xarray[::20]
  167.         y_m = yarray[::20]
  168.         z_m = zarray[::20]
  169.  
  170.     for j in range(0,6):
  171.         print x_m[j], y_m[j], z_m[j]
  172.     #    if (((zarray[j]-n_o)>0) and ((zarray[j+1]-n_o)<0)): #or (((zarray[j]-n_o)<0) and ((zarray[j+20]-n_o)>0)):
  173.     #        counter= counter +1
  174.     #        print zarray[j]-n_o,counter
  175.     #        plotthetaphi(xarray[j],yarray[j],zarray[j])
  176.  
  177.             #ax.plot([xarray[j],xarray[j+20]],[yarray[j],yarray[j+20]],[zarray[j],zarray[j+20]],color='#0000FF' ,marker='.')
  178.             #ax.plot([xarray[j],0],[yarray[j],0],[zarray[j],zarray[j]-n_o],color='#0000FF',marker='.')
  179.             #matplotlib.pyplot.show()
  180.  
  181.  
  182.  
  183.  
  184. def setboundary():
  185.     ax.set_xlim3d(-1,1)
  186.     ax.set_ylim3d(-1,1)
  187.     ax.set_zlim3d(0,1)
  188.  
  189.  
  190. def drawdynamicplot(m,n,d_n, n_o):
  191.     for i in range(0,m):        
  192.             n=vector.rungekutta1(n, d_n)
  193.             d_n=vector.rungekutta2(n, d_n, i)
  194.             setboundary()
  195.             if i%20 == 0:
  196.                 print "-------------" ,i , "----------------"
  197.                 print n[0]                
  198.                 plotpos(n[0],n[1],n[2], n_o, i)
  199.                 plotvel(d_n[0],d_n[1],d_n[2],n[0],n[1],n[2])
  200.                 #if i%320 ==0:               
  201.                 #plotthetaphi(n[0],n[1],n[2])
  202.  
  203. def realtimeconnectline(m,n,d_n):
  204.     for i in range(0,m):
  205.         n=vector.rungekutta1(n, d_n)
  206.         d_n=vector.rungekutta2(n, d_n, i)
  207.         x1=n[0]
  208.         y1=n[1]
  209.         z1=n[2]
  210.         print x1,y1,z1
  211.         if i%10 == 0:
  212.                 x1=n[0]
  213.                 y1=n[1]
  214.                 z1=n[2]
  215.                 if i%20==0:
  216.                     x2=n[0]
  217.                     y2=n[1]
  218.                     z2=n[2]
  219.                 ax.plot([x1,x2],[y1,y2],[z1,z2],color = 'blue',marker='.')
  220.                 matplotlib.pyplot.show()
  221.  
  222.  
  223.  
  224.  
  225. #fig = plt.figure()
  226. #ax = fig.add_subplot(111, projection='3d')
  227.  
  228. c1=-1.0
  229. c2=6.6
  230. c3=-5.4
  231. c4=-c3
  232. c5=1
  233. #[c1,c2,c3,c4,c5]=[-1.0,-2.4,-3.67,-c3,0]
  234. h=.005
  235. vector=vec(1,2,3)
  236. allone=vec(1.0,1.0,1.0)
  237. n=(0.0,.5,.5)
  238. d_n=(5.0,5.0,0.0)
  239. ion()
  240. n=vector.unit(n)
  241. #d_n = vec.unit(vec, d_n)
  242.  
  243. m=5000
  244. xarray=[]
  245. yarray=[]
  246. zarray=[]
  247.  
  248. drawstaticplot(m,n, d_n, n[2])
  249. #drawdynamicplot(m,n, d_n, n[2])
  250. #realtimeconnectline(m,n,d_n)
  251.  
  252.  
Runtime error #stdin #stdout 0.09s 10864KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/8Uccg
language:
Python (cpython 2.7.16)
created:
11 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!