SHARE
TWEET

khí_ly_tưởng

anhhoang2703 Apr 21st, 2017 51 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clc;clear all;close all;
  2. % syms v
  3. kb=1.38e-23;m=6.7e-27; T=300;
  4. v=sqrt((2*kb*T)/m);
  5. L=0.1;
  6. n=1000;
  7. % fv=(m/(2*pi*kb*T))^(3/2)*4*pi*v^2*exp((-m*v^2)/(2*kb*T));
  8. % ezplot(fv,[0 7000])
  9. % shg
  10.  
  11. x=-L+2*L*rand(1,n);
  12. y=-L+2*L*rand(1,n);
  13. goc=-pi+2*pi*rand(1,n);
  14. vx=v*cos(goc);vy=v*sin(goc);
  15. sdeltav=0;
  16. sdeltat=0;
  17. for i=1:1000
  18.  
  19. [Tvc hat trucvc]=mintime(x,vx,y,vy,L);
  20. kvc=length(hat);
  21. % if kvc==1
  22.    
  23.   x=x+vx*Tvc;
  24.   y=y+vy*Tvc;
  25.   if trucvc==1
  26.      vx(hat)=-vx(hat);
  27.      sdeltav=sdeltav+2*abs(vx(hat));
  28.   end
  29.   if trucvc==2
  30.      vy(hat)=-vy(hat);
  31.       sdeltav=sdeltav+2*abs(vy(hat));
  32.   end
  33. % end
  34. % if kvc==2
  35. %      hat;
  36. %      x=x+vx*Tvc;
  37. %      y=y+vy*Tvc;
  38. %      vtg=vx(hat(2));
  39. %      vx(hat(2))=vx(hat(1));
  40. %      vx(hat(1))=vtg;
  41. %      vtg=vy(hat(2));
  42. %      vy(hat(2))=vy(hat(1));
  43. %      vy(hat(1))=vtg;
  44. % end
  45. sdeltat=sdeltat+Tvc;t(i)=sdeltat;
  46. %   plot(x,y,'.',x(2),y(2),'.r','markersize',30);
  47. %
  48. %     axis([-L L -L L]);
  49. %   pause(0.1);
  50. p(i)=(m*sdeltav)/(sdeltat*8*L);
  51.  
  52. end
  53. plt=(n*kb*T)/(4*L^2)
  54. p(end)
  55. plot(t,p,[t(1),t(end)],plt*[1,1]);  
  56.  axis([min(t) max(t) min(p) max(p)])
RAW Paste Data
Top