quantum simulation (collision)

Quantum

    A. First Edition

    This is purely for fun. You see, whenever there is an exam and I feel nervous the best way for me to do is just

    write some free-lance program such that my brain is fully occupied and forget anything at all.

    B.The problem

    The starts long time ago. Once I noticed in Dr. Grogono's web page that these is fancy "snooker" game. And in my

    boring "chem 205" course, the professor mentioned so much about collision of molecular. So, they become inspiration

    for me to write this simple program to simulate the collision of quantum, assuming both moment and energy are

    conserved. In other words, it is the ideal collision.

    To implement quantum, I went such a long way to override various operators of Vector class. The only thing worth

    mentioning is the "*" operator between two vectors. It is weird and I just do it by instinct to reduce my code.

    The physical meaning is not very justified: V(x1,y1)*U(x2,y2) = x1*x2+y1*y2 ==>a number!

    ¡¡

C.The idea of program

    D.The major functions

    The only complicated function is the "collide" which first check if there is possibility of collision by calling

    function of "collisionTime" which returns the time to be collided, -1 indicating impossible to collide.

    Then two colliding quantum "flyto" the position of colliding site where they exchange energy, moment. After that

    collision, both speeds are changed while energy and moment conserves.

    The collisionTime function tries to find the time when distance between two quantum is the sum of their radius.

    (Here I made an interesting mistake. At first, among the two square roots, I choose the positive one and the

    result confused me for quite a long time since the two quantum collides from impossible angle. Then I realized

    that I should choose the smaller root. However, does it matter if I change the original beginning coordinates?)

    I only made a simple test since I am starving now after eight hours' struggle.

    E.Further improvement

    ¡¡

    F.File listing

    1. vector.h

    2. vector.cpp

    3. quantum.h

    4. quantum.cpp

    5. main.cpp

    ¡¡

    file name: vector.h

    class Vector
    {
    protected:
    	double x;
    	double y;
    public:
    	void set(double theX, double theY){x=theX; y=theY;}
    	Vector operator -(Vector& other);
    	Vector& decrement(Vector& other);
    	Vector operator +(Vector& other);
    	Vector& increment(Vector& other);
    	Vector& operator =(Vector& other);
    	Vector& multiply(double scalor);
    	bool operator ==(Vector& other);
    	Vector(double theX=0, double theY=0);
    	double length();
    	double sqrLength();
    	Vector operator *(double scalor);
    	double operator *(Vector& other);
    	void display();
    };

    file name: vector.cpp

    #include <stdlib.h>
    #include <stdio.h>
    #include <math.h>
    #include "quantum.h"

    Vector& Vector::decrement(Vector& other)
    {
    	x-=other.x;
    	y-=other.y;
    	return *this;
    }

    Vector& Vector::increment(Vector& other)
    {
    	x+=other.x;
    	y+=other.y;
    	return *this;
    }

    Vector& Vector::multiply(double scalor)
    {
    	x*=scalor;
    	y*=scalor;
    	return *this;
    }

    Vector::Vector(double theX, double theY)
    {
    	x=theX;
    	y=theY;
    }

    Vector Vector::operator +(Vector& other)
    {
    	Vector result;
    	result.x= x+other.x;
    	result.y= y+other.y;
    	return result;
    }

    Vector Vector::operator -(Vector& other)
    {
    	Vector result;
    	result.x= x-other.x;
    	result.y= y-other.y;
    	return result;
    }

    Vector& Vector::operator =(Vector& other)
    {
    	x= other.x;
    	y= other.y;
    	return *this;
    }

    double Vector::operator *(Vector& other)
    {
    	return x*other.x+y*other.y;
    }

    Vector Vector::operator *(double scalor)
    {
    	Vector result;
    	result.x = x*scalor;
    	result.y = y*scalor;
    	return result;
    }

    bool Vector::operator ==(Vector& other)
    {
    	return x==other.x && y==other.y;
    }

    double Vector::sqrLength()
    {
    	return (x*x+y*y);
    }

    double Vector::length()
    {
    	return sqrt(x*x+y*y);
    }

    void Vector::display()
    {
    	printf("(%f, %f)", x, y);
    }

    file name: quantum.h

    #include "vector.h"

    const double DefaultStartMoment=0;
    const double DefaultRadius=1;
    const double DefaultSpeed=5;
    const double DefaultStartPos=0;

    //assume the quantum always has same mass, it is much simple
    class Quantum
    {
    protected:
    	Vector coord;
    	Vector speed;
    	double elapsed;
    	double radius;
    	void initialize();
    public:
    	Quantum(Vector& position, Vector& theSpeed, double theMoment, double theRadius);
    	Quantum();
    	void setRadius(double newRadius){radius=newRadius;}
    	void setSpeed(double vx, double vy){speed.set(vx, vy);}
    	void setSpeed(Vector& speedVector){speed=speedVector;}
    	void setPos(double px, double py){coord.set(px, py);}
    	void setPos(Vector& posVector){ coord=posVector;}
    	void setMoment(double now){elapsed=now;}
    	double collisionTime(Quantum& other);//return the time
    	Vector fly(double elapse);//the pos vector quantum will fly to
    	void flyto(double elapse);
    	void collide(Quantum& other);
    	void display();
    };

    ¡¡

    file name: quantum.cpp

    #include "quantum.h"
    #include <stdlib.h>
    #include <stdio.h>
    #include <math.h>

    Quantum::Quantum()
    {
    	initialize();
    }

    void Quantum::initialize()
    {
    	radius=DefaultRadius;
    	elapsed=DefaultStartMoment;
    	coord.set(DefaultStartPos, DefaultStartPos);
    	speed.set(DefaultSpeed, DefaultSpeed);
    }

    Quantum::Quantum(Vector& thePos, Vector& theSpeed, double theMoment, double theRadius)
    {
    	coord=thePos;
    	speed=theSpeed;
    	elapsed=theMoment;
    	radius=theRadius;
    }

    Vector Quantum::fly(double elapse)
    {
    	return coord+(speed*elapse);
    }

    /*
    assume P{(x1,y1),(vx1,vy1)}, Q{(x2,y2),(vx2,vy2)}
    flying time t
    P{(x1+vx1*t, y1+vy1*t)}, Q{(x2+vx2*t),(y2+vy2*t)}
    The distance between P,Q should be (r1+r2)
    So, square(x2-x1+t(vx2-vx1))+square(y2-y1+t(vy2-vy1))=square(r1+r2)
    a=square(vx2-vx1)+square(vy2-vy1)
    b=2*((vy2-vy1)(y2-y1)+(vx2-vx1)(x2-x1))
    c=square(x2-x1)+square(y2-y1)-square(r1+r2)
    */

    double Quantum::collisionTime(Quantum& other)
    {
    	double a,b,c, root;
    	a=(speed-other.speed).sqrLength();
    	//here is the * of two vector, like angle between them
    	b=((speed-other.speed)*(coord-other.coord))*2;
    	c=(coord-other.coord).sqrLength()-(radius+other.radius)*(radius+other.radius);
    	root=b*b - 4*a*c;
    	if (root<0)
    	{
    		return -1;//won't collide at all
    	}
    	return -(sqrt(root)+b)/(2*a);
    }

    void Quantum::flyto(double elapse)
    {
    	coord.increment(speed*elapse);
    	elapsed+=elapse;
    }

    /*
    assume the collision point is p1(x1,y1), p2(x2,y2)
    the delta speed vector should be on the line of p1_p2
    delta=scalor*(x1-x2,y1-y2) is exerted on Quantum Q
    -delta should be exerted on Quantum P
    the conservation of moment is satisfied since |delta|=|-delta|
    or delta+(-delta)=0 moment
    the conservation of energy should also be satisfied
    sq(vx1)+sq(vy1)+sq(vx2)+sq(vy2)=sq(vx1+scalor*(x1-x2))+sq(vy1+scalor*(y1-y2))
    +sq(vx2+(scalor*(x2-x1))+sq(vy2+(scalor*(y2-y1))
    simplify and get:
    a=(sq(x1-x2)+sq(y1-y2))
    b=(x1-x2)(vx2-vx1)+(y1-y2)(vy2-vy1)
    scalor=b/a
    */

    void Quantum::collide(Quantum& other)
    {
    	double elapse, scalor;
    	Vector delta;
    	if ((elapse=collisionTime(other))<0)
    	{
    		return;
    	}
    	flyto(elapse);
    	other.flyto(elapse);
    	//these are speed delta
    	delta=coord-other.coord;
    	scalor= ((other.speed-speed)*delta)/(delta.sqrLength());
    	speed.increment(delta*scalor);
    	other.speed.decrement(delta*scalor);
    }

    void Quantum::display()
    {
    	printf("coord:");
    	coord.display();
    	printf(" speed:");
    	speed.display();
    	printf("\n");
    }

    ¡¡

    file name: main.cpp

    #include "quantum.h"
    #include <stdlib.h>
    #include <stdio.h>

    int main()
    {
    	Quantum Q, P;
    	//double moment;
    	Q.setPos(1,2);
    	Q.setSpeed(1,0);
    	P.setPos(4,1);
    	P.setSpeed(-2,1);
    	printf("before collision:\n");
    	printf("now Q is:");
    	Q.display();
    	printf("now P is:");
    	P.display();
    	/*
    	if ((moment=Q.collisionTime(P))>0)
    	{
    		printf("Q and P will collide at %f\n", moment);
    	}
    	P.flyto(moment);
    	Q.flyto(moment);
    	*/
    	Q.collide(P);

    	printf("after collision:\n");
    	printf("now Q is:");
    	Q.display();
    	printf("now P is:");
    	P.display();
    	return 0;
    }

    ¡¡

running result:

before collision:
now Q is:coord:(1.000000, 2.000000) speed:(1.000000, 0.000000)
now P is:coord:(4.000000, 1.000000) speed:(-2.000000, 1.000000)
after collision:
now Q is:coord:(1.367544, 2.000000) speed:(-2.000000, 1.000000)
now P is:coord:(3.264911, 1.367544) speed:(1.000000, 0.000000)
Press any key to continue

¡¡