气体分子小球模型Processing代码运行界面
参考文章：Experiments in statistical mechanics
Main function
float ball_l=80;
int N = 3;
Particle [] ps= new Particle[N];
Bar bar = new Bar (int(ball_l)*5);
float heigh = ball_l*4;
float widt = ball_l*8;
//parameters for statistics
float timeBoxWidth=ball_l*4; //the width of time box
float timeBoxHeight=ball_l*2;
float timeBoxOriginX=ball_l*8; //leftupper corner x of time box
float timeBoxOriginY=ball_l*0; //leftupper corner y of time box
int counter = 0; //time step counter
float [] ts = {}; //time steps
float [] ys = {}; //bar locations
float histBoxWidth = ball_l*4; //the width of histogram box
float histBoxHeight = ball_l*2; //the height of histogram box
float histBoxOriginX=ball_l*8; //leftupper corner x of hist box
float histBoxOriginY=ball_l*2; //leftupper corner y of hist box
void setup() {
size(int(ball_l)*12,int(ball_l)*4);
/*
for (int i=0; i<ps.length; i++){
ps[i] = new Particle(new PVector(random(width),random(height)),ball_l/2);
}
*/
ps[0] = new Particle(new PVector(ball_l*1.5,ball_l*1.5),ball_l/2);
ps[1] = new Particle(new PVector(ball_l*6.5,ball_l*3.5),ball_l/2);
ps[2] = new Particle(new PVector(ball_l*7.5,ball_l*0.5),ball_l/2);
//ps[3] = new Particle(new PVector(ball_l*7.5,ball_l*1.5),ball_l/2);
}
void draw() {
background(240);
counter+=1;
//plot lattice
for (int x=0; x<=8; x++){
line(x*ball_l,0,x*ball_l,heigh);
}
for (int y=0; y<4; y++){
line(0,y*ball_l,widt,y*ball_l);
}
//display particles
for (Particle p: ps){
p.run();
}
//display bar
bar.run();
//collision between particles
for (int i=0; i<ps.length; i++){
for (int j=i+1; j<N; j++){
PVector dl = PVector.sub(ps[j].location, ps[i].location);
float dis = dl.mag();
if (dis < (ps[i].r+ps[j].r)){
//balls have contact so push back
//https://channel9.msdn.com/Series/Sketchbooktutorial/Simple-Collision-Detection-and-Response
PVector normal = PVector.div(dl,dis);
PVector midpoint = PVector.div(PVector.add(ps[i].location, ps[j].location),2);
ps[i].location = PVector.sub(midpoint,PVector.mult(normal,ps[i].r));
ps[j].location = PVector.add(midpoint,PVector.mult(normal,ps[j].r));
float scale = PVector.sub(ps[i].velocity, ps[j].velocity).dot(normal);
PVector dv = PVector.mult(normal,scale);
ps[i].velocity.sub(dv);
ps[j].velocity.add(dv);
}
}
}
// collision between particles and the bar
for (int i=0; i<ps.length; i++){
float dx = bar.x - ps[i].location.x;
if (abs(dx) < ps[i].r){
//print (dx+"___");
//contact so push back
float normalX;
if (dx > 0){normalX = 1;}
else {normalX=-1;}
float midpointX = (ps[i].location.x+ps[i].r*normalX+bar.x)/2;
ps[i].location.x = midpointX - normalX*ps[i].r;
bar.x = midpointX;
float dvx = ps[i].velocity.x - bar.velocity;
ps[i].velocity.x -= dvx;
bar.velocity += dvx;
}
}
//calculate total energy
float E=0;
for (int i=0; i<ps.length; i++){
E+=sq(ps[i].velocity.mag());
}
E+=sq(bar.velocity);
//print (E+"_____");
//-------------------time box-------------------
fill(255);
rect(timeBoxOriginX,timeBoxOriginY,timeBoxWidth,timeBoxHeight);
fill(100);
text("Time = "+counter,timeBoxOriginX+timeBoxWidth-80,timeBoxOriginY+timeBoxHeight-10);//print counter
text("Bar location",timeBoxOriginX+10,timeBoxOriginY+20);
//add a new data point
ts=(float[]) append(ts,counter);
ys=(float[]) append(ys,bar.x);
//obtain rescaled xs and ts
float [] rts = new float[ts.length];
float [] rys = new float[ys.length];
for (int i=0; i<ts.length; i++){
rts[i]=timeBoxOriginX+timeBoxWidth*ts[i]/( max(ts)+1);
rys[i]=timeBoxOriginY+timeBoxHeight*(1-ys[i]/( max(ys)+1));// change y direction
}
//draw lines connecting time points
if (rts.length>1){
for (int i=0;i<rts.length-1; i++){
float x1=rts[i];
float x2=rts[i+1];
float y1=rys[i];
float y2=rys[i+1];
strokeWeight(0.5);
fill(100,100);
line(x1,y1,x2,y2);
//point(x1,y1);
fill(200);
}
}
//-------------------histogram-------------------
fill(255);
rect(histBoxOriginX,histBoxOriginY,histBoxWidth,histBoxHeight);
fill(100);
text("Bar location distribution",histBoxOriginX+10,histBoxOriginY+20);
int Nbars = 20;
int [] histData = new int[Nbars];
float bin = histBoxWidth/Nbars;
for (int i=0;i<ys.length; i++){
int n = int(Nbars*(ys[i]-min(ys))/( (max(ys)-min(ys))+0.001)); //the ith element belongs to the nth bin
histData[n]+=1;
}
for (int i=0;i<Nbars;i++){
float h = histBoxHeight*histData[i]/float(ys.length);
rect(histBoxOriginX+i*bin,histBoxOriginY+histBoxHeight-h,bin,h);
}
Particle Class
class Particle {
PVector location;
PVector velocity;
//PVector acceleration;
float mass;
float r;
Particle(PVector l, float radius) {
mass = 1;
//acceleration = new PVector(0,0);
velocity = new PVector(5,5);
location = l.get();
r = radius;
}
void run() {
update();
display();
checkEdges();
}
void update() {
//velocity.add(acceleration);
location.add(velocity);
}
void display() {
stroke(0);
fill(175);
ellipse(location.x,location.y,r*2,r*2);
}
void checkEdges() {
if (location.x < r){
location.x = r;
velocity.x *= -1;
}
if (location.x > widt-r){
location.x = widt-r;
velocity.x *= -1;
}
if (location.y < r){
location.y = r;
velocity.y *= -1;
}
if (location.y > heigh-r){
location.y = heigh-r;
velocity.y *= -1;
}
}
}
Bar class
class Bar {
float x;
float velocity;
//float acceleration;
float mass;
Bar(float l_) {
mass = 1;
//acceleration = 0;
velocity = 0;
x = l_;
}
void run() {
update();
display();
checkEdges();
}
void update() {
//velocity+=acceleration;
x+=velocity;
}
void display() {
stroke(0);
fill(100);
strokeWeight(2);
line(x,0,x,heigh);
strokeWeight(1);
}
void checkEdges() {
if (x < 0){
x = 0;
velocity *= -1;
}
if (x > widt){
x = widt;
velocity *= -1;
}
}
}
}
网友评论