Section 1: A simple Vector class
Part A: Creating the class
Classes are a way in C++ to combine data and functions to make a new user-defined “type”. Examples of types that you have already used that are built using classes are things like if stream, of stream and string. Also, you have used some user defined types that aren’t in the standard library when doing the exercises from SectionC – for example the SPA::Window class.
Your first step in this worksheet should be to go and read the document from Section D on classes – “SimpleClassIntroduction”. Then you’ll start writing your own class to implement a simple 3D vector.
Once you have read the document carry on with this worksheet.
In the first part of the worksheet you’ll create a very basic class for storing and manipulating three values that represent Cartesian coordinates of a 3D vector. In Clion there are 2 ways to create the files for a new class. The simplest is to right-click on your project and select New->C++ Class.
In the window that appears enter the name of the class at the top (no spaces allowed) and be sure that the tick-box labelled “Add to targets” is selected and that the first box of the two further checkboxes is ticked. The other way is to add a cpp file in the same way that you did for functions.cpp in SectionC but then you’ll have to type the class declaration in yourself. The former method will insert this code for you.
- Create a new project called Ex3.2
- Create a new class called Vector
Note: be careful with the capitalisation of your class, there is already a vector class that you have encountered called std::vector. If you call your class vector instead of Vector you could get a naming clash. There is a feature of C++ that we have already encountered to deal with these kind of naming clashes – an extra mark will be granted if you can work out what it is and use it correctly in this project.
- In your class add three private data members. They should all be doubles and be called x, y and z. These will represent the Cartesian coordinates of a 3D vector. Using the examples in “SimpleClassIntroduction” (or from elsewhere) add 6 functions to your class – one to set the value and one to get the value of each of the three coordinates. For example: void setX(double v) and double getX(). Hint: Remember functions that don’t modify any of the data of the class should be declared “const”. The const should be used in both the declaration and definition (if the definition is in a
- Using the examples in “SimpleClassIntroduction” (or from elsewhere) add 6 functions to your class – one to set the value and one to get the value of each of the three coordinates. For example: void setX(double v) and double getX(). Hint: Remember functions that don’t modify any of the data of the class should be declared “const”. The const should be used in both the declaration and definition (if the definition is in a
- Hint: Remember functions that don’t modify any of the data of the class should be declared “const”. The const should be used in both the declaration and definition (if the definition is in a CPP file )
- Add a constructor to your class that takes three arguments and then using the initializer syntax described in the “SimpleClassIntroduction” document, add initializers for the 3 variables belonging to your class. Add default values that set the components to 0 if no arguments are given.
- Commit and push your code now.
Part B: Testing your code
A critical part of programming is testing your code as you go along. If you can write your code in relatively self-contained components and test each one as you go along it is far easier to debug and be sure that your code is doing the correct thing.
When comparing real numbers (ie those with decimal places) in code it is always important to realise that there is no perfect representation of every possible number and so although the results of two different calculations might be the same mathematically, depending on the exact calculations they might not be exactly the same in computer code. This can make it difficult to test if two floating point numbers are actually the same. To make this easier you will now write a little helper function.
- Create a new source code file, functions.cpp and the associated header file functions.h
- In functions.h write a function declaration for the following function:
bool testFloat(double value,double test,double p=0.0000001)
- In functions.CPP write the appropriate function definition. The function body should carry out the following steps:
oCalculate the difference between value and test and store this in a double called diff
o Using an if statement test that the absolute value (ie ignoring the sign) of diff is less than p
o Return a value of true if this is the case
o Return a value of false if it is not
- Hint: remember a default value should only appear in the declaration – you wil get an error if it appears in the definition as well.
Now you will write a test function that makes use of this together with your new Vector class.
- In functions.h write a function declaration for the following function:
In functions.cpp write the appropriate function definition with a function body that carries out the following actions:
oCreate a new Vector object
o Test (using testFloat with p=0.0001) that the initial x,y and z values are 0.
o Using the appropriate member function set the value of x to 10.
o Using the appropriate member function set the value of y to 20.
o Using the appropriate member function set the value of z to 30.
o Test (using testFloat with p=0.0001) that the x y and z values of your Vector object are correctly set to 10., 20., and 30.0 respectively.
o Your function should return true if all the tests are correct and false if any of them fail.
Now it’s time to put this together to do a full test of your code
- In main.cpp write code to carry out the following steps:
o Call the function test1() and store the return value in a bool called t1
o Using an if statement test the value of t1 and carry out the following actions if it is true:
o Print the string “Test 1 SUCCESS” to the console
o else carry out the following steps if it is false:
o print the string “Test 1 FAILED” to the console
Now compile and test all your code. Make sure that your test code is correct and then run it. If you don’t see “Test 1 SUCCESS” go back and see if you can figure out where the code is going wrong. Adding extra cout statements to your code is a good way to check what values are being calculated at different steps of the program.
- Commit and push your code now
Part C: Adding more methods
Currently the vector class is kind of bare, it just stores the data and provides ways to retrieve and change it. To make the Vector class more useful you need to add some more functions to perform some useful operations.
Add the following methods to your Vector class – put the declarations in the h file and the definitions in the cpp file
- modulusSqr : add a const member function with this name. It should return a double equal to the squared modulus of the vector:
|?|² = ? ² + ? ² + ? ²
- modulus: add a const member function with this name. It should return the modulus of the vector.
- dot : add a const member function with this name. It should take on argument that should be a const reference to a Vector object. This function should return a double that is equal to the dot product of the Vector it is being called on and the vector in the argument.
- Hint: to access the data members of an object inside a member function you just use their names. To access the data members of an argument to the function you use the usual access methods:
o Remember the dot product of two vectors u and v is defined as:
o distance : add a const member function with this name. It should take an argument that should be a const reference to a Vector object. This function should return a double that is equal to the distance between two Vectors – assuming they are position vectors. The square of the distance is defined as:
o Hint: create a new Vector object initialised to the difference between the two vectors in vector form and return its modulus.
Commit and push your code now!
Now you need to add some test code to make sure these new methods are working correctly.
Add new functions to functions.h/function.cpp (in each case make use of testFloat where appropriate)
o testModulus: add a function with this name that returns a bool and takes no arguments. In the body of the function write code to create a vector (1,2,3).
- Test that the modulus and modulusSqr functions return the correct value using testFloat with p=0.0001
- The function should return true if both modulus and modulusSqr return the correct value
o Add code to main.cpp to call this new function and print “Test Modulus : SUCCESS” if it returns true and “Test Modulus : FAILED” if it returns false.
o testDistance: add a function with this name that returns a bool and takes no arguments. In the body of the function write code to create a vector v1(1,2,3) and a vector v2(3,2,1). Test that a call to the distance function returns the correct value using testFloat with p=0.0001.
o Hint: getting the distance should look something like:
double d = v1.distance(v2);
o Add code to main.cpp to call this new function and print “Test Distance : SUCCESS” if it returns true and “Test Distance : FAILED” if it returns false.
otestDot: add a function with this name that returns a bool and takes no arguments. In the body of the function write code to create a vector v1(1,0,1) and a vector v2(0,1,0). Test that your dot member function correctly calculates the dot product of these two vectors – it should be 0. Make use of testFloat with p=0.0001
oAdd a further test – use your dot function to calculate the dot product of v3(1,1,1) and v4(0,1,1) it should be equal to 2.
o Your function should return true if both these tests are correct and false if either or both of them fail.
o Add code to main.cpp to call your new function and print “Test Dot : SUCCESS” if it returns true and “Test Dot : FAILED” if it returns false.
Commit and push your code.
Check your test functions very carefully! If they are incorrect you might think your other code is working correctly when it is not. Make sure all the tests work and print SUCCESS.
Part D : Some more complicated methods!
Now there are three more methods needed before moving onto class operators:
scale: add a function with this name to your class. It should have 1 double argument f and return nothing. The body of the function should multiply each of the data members x, y and z by the value of f. Note since this function modifies the values of the data it cannot be labelled as const.
unit: add a function with this name to your class. It should take no arguments and return a Vector object. In the body of the function create a new Vector v with x, y and z arguments equal to the data members x, y and z. Then scale this vector using your scale function such that it has unit modulus. The function should then return this scaled vector as its return value.
cross: add a function with this name to your class. It should take one const Vector reference argument and return a Vector object. In the body of the function calculate the cross product of the object on which the function is called and the vector argument. Return the resulting vector as the return value. Remember the cross product ? of two vectors ? and ? is:
Now add some more test code for each of these new methods – each test function should, as before, take no arguments and return a bool. In each case use testFloat where appropriate with p=0.0001
o testScale : add this function to functions.cpp/h. In the body of the function create a new vector v(1,1,1) and scale it with the scale method by a factor of 3. Test that the x, y and z values of the scaled vector are each equal to three. The function should return true if this is the case and false otherwise.
o Add a call to this function in main.cpp following the same patterns as the previous tests. Print “Test Scale : SUCCESS” if your function returns true and “Test Scale : FAILED” otherwise.
o testUnit : add this function to functions.cpp/h. In the body of the function create a new vector v1(3,3,3). Assign the return value of calling the unit method to another vector called v2. Test that the modulus of v2 is equal to one. Then test that the dot production of the two vectors is equal to the modulus of v1. Your function should return true if both these cases are correct and false otherwise.
o Add a call to this function in main.cpp following the same pattern as the previous tests. Print “Test Unit : SUCCESS” if your function returns true and “Test Unit : FAILED” otherwise.
o testCross : add this function to functions.cpp/h. In the body of the function:
- create a vector v1(1,0,0), a vector v2(0,1,0) and a vector v3(0,0,1).
- Calculate v1.cross(v2) and assign the result to a new vector called v4.
- Calculate v2.cross(v3) and assign the result to a new vector called v5.
- Calculate v3.cross(v1) and assign the result to a new vector called v6.
- Test that v4 = (0,0,1)
- Test that v5 = (1,0,0)
- Test that v6 = (0,1,0)
- Hint: use the distance function and the appropriate vectors v1,v2 and v3
- Return true if all the tests are correct, otherwise return false
o Add a function call to main.cpp to call this new function. Print “Test Cross : SUCCESS” if it returns true and “Test Cross : FAILED” if it returns false.
You now have an almost fully functioning vector class. Commit and push your work!
Part E : Final steps in defining the class – adding the operators
Before you proceed take a look at the document: “MoreClasses”. Amongst other things this document discusses how to define class operators.
In the part you will add the basic operators for adding, subtracting and inverting vectors.
The most compact way to do this is first to define the composite operators -= and +=:
Remember: a += b; is equivalent to: a = a +b;
Here is an example on how to do this for class with a single data member:
Add this method to your class – you can put the definition in the .h file in this case.
Note that it is not const since it is modifying the data members of the class. Extend this to operate correctly for all three data members.
Now you can add the operator+ that defines the operation of addition in terms of this composite operator. It should look like this:
Here you are creating a new vector v2 – that is a copy of the vector it is being called on – adding the vector in the argument to it using the composite addition assignment operator and returning the resulting value. Remember from the MoreClasses document code that looks like: a + b is equivalent to a.operator+(b) .
- Put this example into your class and also implement the equivalent operators for operator-= and operator-.
There is one final operator to add and that is “unary negation”. This is basically the unary minus operator, you need it to be able write code like:
Vector v2 = -v1;
Implement the unary negation operator in your class. It should take no arguments and return a Vector. In the body of the function you should create a new vector whose x, y and z values are the negative of the object’s x, y and z values. The method should be const.
Now for some test code
Add a new function to functions.cpp/h called: testOperators. As with the other test functions it should take no arguments and return a bool. It should contain three basic tests:
o Test 1:
- Create a vector v1(1,2,3)
- Create a vector v2(2,4,6)
- Create a vector v3=v1+v1
- Test that v3 is equal to v2
- Hint: test its distance from v2 is 0 using testFloat with p =0.0001
o Test 2:
- Create a vector v4 = v2-v1
- Test that v4 is equal to v1
o Test 3:
- Create a vector v5 = -v1
- Create a vector v6(-1,-2,-3)
- Test that v5 = v6
o Your function should return true if all the tests succeed, otherwise it should return false
Add a call to your new function in main.cpp and print “Test Operators : SUCCESS” if it returns true and “Test Operators : FALSE” if it returns false
Section 2: Adding a Particle class
In this section you’ll add a particle class with a collection of useful methods to build a simulation of colliding particles in a box.
Part A: Basic particle class
- Create a new class in your project called Particle
- Add two data members of type Vector called p and v to represent the position and velocity of the particle
- Add one data member of type double called radius to represent the radius of the spherical particle.
- Add a constructor that takes three arguments:
- A double called r
- A const reference to a vector called p – the position of the particle
- A const reference to a vector called v – the velocity of the particle
- Use the correct form of the initializer syntax from “MoreClasses” to set the radius, position and velocity of the particle correctly.
- Use default arguments in the constructor such that if the constructor is called without arguments the position is (0,0,0) and the velocity is (0,0,0) and radius is 1.
- Hint: You can create a default Vector with Vector()
- Add four methods to your class to allow you to set the position and velocity and to retrieve the position and the velocity. These should be called: setPosition, setVelocity and getPosition and getVelocity. Use the appropriate return types and const where required
- Add a new test function to functions.cpp/h called testParticle. This should take no arguments a return a bool.
- Add a further method to retrieve the radius called getRadius. This should take no arguments, be const and return a double equal to the value of radius.
- In this function write code to create a new particle with position (0,0,0) and velocity (0,0,0) using the default contructor (ie constructor with no arguments)
- Using your getPosition and getVelocity test that the particle has been initialized correctly.
- Using setPosition and setVelocity set the position and velocity of the particle to (1,2,3) and (4,5,6) respectively. Again using getPosition and getVelocity test that the particle has the correct position and velocity.
- In all your tests make appropriate use of testFloat with p =0.0001
- Your function should return true if all tests pass and false if they fail.
- Add the appropriate line to main.cpp to call your function, printing “Test Particle : SUCCESS” if it returns true and “Test Particle : FAILED” otherwise.
Part B : Moving and Colliding Particles
In this part you will continue working on your particle class – adding code to make the particles move and collide.
Add a new method to your Particle class:
o move : the member function should have this name. It should return a void and take one argument – a double called t. It is going to change the position of the particle so it cannot be const. Put the definition of the function in Particle.cpp. The body of the function should move the particles position assuming t is the time in appropriate units. You can implement this in a way that is independent of the Cartesian nature of the vector representation we are using in the following way:
- create a new Vector called delta as a copy of v – the particle velocity.
- scale delta by t
- add delta to the position vector member of the particle
o This approach hides the vector representation in the vector class
Add a new test function to functions.cpp/h:
o testMove : the function should have this name, take no arguments and return a bool.
o In your test create a new particle with position (0,0,0) and velocity (1,1,1).
o Then call the move function with time = 1
o Test that the new position of the particle is (1,1,1)
o Your function should return true if the test succeeds and false otherwise.
o Add the appropriate code following the established pattern to call your test function in main.cpp and print “Test Move : SUCCESS” if your function returns true and “Test Move : FAILED” otherwise
The next step is to add code to work out the time to collision between two particles:
The easiest way to consider when two particles will collide is to consider their relative positions and velocities:
From this you know the variation of r with time, t:
You can calculate the time of closest approach in two ways – either by taking the derivative of the modulus squared of r(t) with respect to t and finding where it reaches 0, or more simply, consider
that at the point of closest approach the relative positions of the particles is perpendicular to the velocity:
And thus the time at distance of closest approach is:
If this is negative then, the particles are heading away from each other and so won’t collide.
This isn’t quite the time to collision though – if the distance of closest approach is too big the particles won’t meet.
To calculate the distance of closest approach, b, use:
b²=r(0)² – s²
Where s is just the distance travelled in time ??? :
Now, if the distance of closest approach is less than twice the radius of the particles, they will collide – otherwise they will pass by each other.
If the distance of closest approach is exactly twice the radius the the time of closest approach is also the time of collision.
At the time of collision you know that the particles must be twice the radius, R, apart – so the distance, d, prior to the point of closest approach is:
d²=4r² – b²
And the time to travel that distance, is just d/|v|, so the overall time to collision is just the time to closest approach minus this time:
Implement a new function in your Particle class called timeToCollision. This function should return a double and take one argument as a const reference to a Particle called particle. The function definition should go in Particle.cpp and implement the calculation described above in the following steps:
o Create a new vector r0 that is the difference between the particle position in the argument and the particle position in the object being called. For example if you were to call the function somewhere using two particles p1 and p2 with position vectors r1 and r2, then r0 = r2-r1.
double t = p1.timeToCollision (p2) ;
oCreate a new vector v that is the relative velocity of the two particles: v2-v1
o Using the equation for ??? above calculate the time of closest approach of the two particles and store it in a double called tca.
o Test if tca is less than zero – if so the particles will not collide in the future – return a value of FLT_MAX (this is defined in <cfloat>) and exit the function.
o Create a new Vector called s that is a copy of the relative velocity and scale using the scale function by ???.
o Evaluate the distance of closest approach squared using the equations above and test if it greater than 4 times the radius squared – if so the particles will never collide return a value of FLT_MAX from the function.
o Otherwise – calculate the quantity ???? described above and return this as the result of the function.
Now you need to test this function by constructing some collision scenarios. Add a test function to functions.cpp/h called testCollision, taking no arguments and returning a bool.
This function should carry out the following tests and only return true if all succeed:
o Test 1: Head on collision
o Create two particles p1 and p2
o p1 : radius = 1, position = (0,0,0), velocity = ( 1,0,0)
o p2: radius = 1 , position = (12,0,0), velocity = (-1,0,0)
o Since the particles have radius 1 a collision should occur after 5 seconds.
o Test 2: Non-head-on collision
o Create two particles p3 and p4
o P3: radius = 1, position = (0,0,0), velocity = (1,1,0)
o P4: radius = 1, position = (5,5,6.99999), velocity = (0,0,-1)
o Again the time to collision should be 5 seconds
o Test 3: No collision
oCreate two particle p5 and p6
oP5: radius = 1, position = (0,3,0), velocity = (1,0,0)
oP6: radius = 1, position = (12,0,0), velocity = (-1, 0 , 1)
oTest that the result is equal to FLT_MAX
Part C : Implementing a particle-particle collision
Now you need to implement the full collision – this includes computing the time to collision, moving the particles and adjusting the velocities.
Add a new member function to your Particle class called collideParticles. It should return nothing and take a single reference argument to a Particle. Note this shouldn’t be a const reference because you will need to be able to modify the particle.
Put the function definition in Particle.cpp implementing the following logic:
o First calculate the time to collision between the particles.
o If the result is less than 0 or very close to FLT_MAX then just return from the function without doing the collision
o Move both particles to the collision point using the move function. Remember this is a member function so the first particle is the one the function belongs to and the second particle is the one passed in the argument.
o When two spheres collide they exchange velocities along the vector joining their centres.
o Get a unit vector along the line joining the two particles (ie along r2-r1)
o Make two copies of the unit vector called v1 and v2.
o Scale v1 by the dot product of the velocity of particle 1 and n to get the velocity vector along n of particle 1.
o Scale v2 by the dot product of the velocity of particle 2 and n to get the velocity vector along n of particle 2.
o Calculate the change in velocity by taking the different v2-v1 and storing it in a new Vector called deltaV
o The new velocity vector for particle 1 is v+deltaV
o The new velocity vector for particle 2 is v-deltaV
o Where v is the original velocity of the particle in each case
Now add a test function for this code:
o Add a new function to functions.cpp/h called testCollision. Again, it should take no arguments and return a bool.
o In the body of the function do the following 2 tests:
o Create a new particle p1: radius=1, position = (0,0,0), velocity = (1,0,0)
o Create a new particle p2: radius=1 position = (5,0,0), velocity = (0,0,0)
o Cause the particles to collide by calling p1.collideParticles(p2)
o Now check that the velocity of p1 is (0,0,0) and the velocity of p2 is (1,0,0)
o For the second test:
o Create a new particle called p3: radius =1, position =- (0,0,0), v(1,0,0)
o Create a new particle called p4:radius=1,position (10,0,0),v(-0.5,0,0)
o Collide the particles and check the x-component of the velocity of p1 is now -0.5 and the x-component of the velocity of p2 is now 1.
Part D : Putting the particles in a box
Now for the final steps: putting the particles in a box. Imagine a cubical box whose centre has coordinates (0,0,0).
Create a new class called Box. Box should have two data members: a std::vector of Particles called particles and double called size. The vector represents the list of particles in the box, and size represents the side length of the cube (so given the centre is 0,0,0 each side is 250 units away from the middle).
Add a constructor that takes one double argument used to initialize the value of size. Use a default argument value of 500.
The class should have three member functions
o addParticle : this should return nothing and take a single argument of a const reference to a Particle
- Implement this method in Box.cpp it should just add the Particle in the argument to the vector of Particles.
o timeToWallCollision: this function will work out the time to the next wall collision for a particle
- implement this method in Box.cpp. It should return a double and take two arguments. The first argument should be an integer representing the wall to be struck. The second argument should be a const reference to a Particle. The body of the function should be implemented in Box.cpp and should do the following:
o If the wall argument is 1 consider the x-coordinate, 2 consider the y-coordinate and 3 consider the z-coordinate. Calculate the time to hit the wall in that particle direction by computing the following (example given for x-component only):
o If r is the position vector of the particle and v is the velocity vector and s is half the side length of the box, then calculate:
o doWallCollision: this function should take two arguments: the first an integer representing the wall to be struck. The second reference to a particle. Implement the function in Box.cpp. It should carry out the following actions:
- calculate the time for the particle passed as argument to hit the wall using the timeToWallCollision function. Move the particle by that amount of time and then invert the component of velocity of the particle corresponding to the wall struck (1 =X, 2=Y, 3=Z). This represents a collision of a particle with a wall of significantly higher mass.
o Step : this function should return a double and take no arguments.
This function will loop over all particles and determine what type of collision will occur next for that particle and when.
Then it should carry out the collision that happens first and return the time taken for that to happen.
Implement this function in Box.cpp as follows:
- Create a new double called timeToNextCollision and initialise it to FLTMAX
- Create a new integer called collisionType and initialize it to 0
- Create a new integer called collider1 and initialize it to -1
- Create a new integer called collider2 and initialize it to -1
- Now loop over all particles in the Particle vector.
- For each particle calculate the time for it to collide with wall 1,2 and 3
- Take whichever is the shortest time and compare it with timeToNextCollision. If it is shorter than timeToNextCollision do the following:
o Store the time in timeToNextCollision
o Set collider1 equal to the index of the current particle in the vector
o Set the collisionType equal to the wall number
This determines which particle will hit a wall next and which wall it will be.
Now loop over all pairs of particles.
Hint: use two nested loops to cover each pairing only once
For each pair calculate their time to collision using the appropriate method of the Particle class.
If the time is less than timeToNextCollision then do the following:
o Set timeToNextCollision to the time
o Set the collisionType to 0
o Set collider1 to the index of the first particle
o Set collider2 to the index of the second particle
The next step is to move forwards in time by “timeToNextCollision”
First carry out what ever collision is appropriate:
- If collisionType is 0, then use Particle::colliderParticles for the particles given by the indices stored in collider1 and collider2
- Otherwise, use Box::doWallCollision for the appropriate wall stored in collisionType and Particle with index collider1.
Then loop over all particles
- If the index of the particle is collider1 or collider2 then do nothing since these particles have already been handled. Otherwise move the particles using Particle::move with the appropriate time parameter.
Finally, return from Step() with the value of the amount of time that elapsed for that collision ie timeToNextCollision.
o Now implement code in main.cpp to do the following:
- Look at the instructions in SectionDBasicWorksheet PartD for instructions on how to get hold of the code to generate normally distributed random numbers.
- Read the worksheet on RandomNumbers
- Create a Box object in your main.cpp
- In a loop create a particle with radius 1, and with a uniform but random spatial distribution within the box (size 500 centred on (0,0,0)) and x, y and z components of velocity chosen from random normally distribution numbers with a mean of 0 and a standard deviation of 10.
- Add the particle to the Box.
- Loop 100 times to put 100 particles in the box.
- Then in a new loop, carry out 10000 collisions by calling the Step function.
o You now have a multiple particle kinetic simulation.
o Some things to investigate if you have time:
What is the total kinetic energy of the particles in the box?
- Does this appear to be conserved by your simulation?
What is the total momentum of the particles in the box?
- Is this conserved? (Hint. It isn’t, but can you explain why?)
How does the “speed” distribution evolve with time?
How might you extract this data to do further analysis?