Skip to content
Snippets Groups Projects
Commit 38fe6fe4 authored by cprutean's avatar cprutean
Browse files

Update 2 files 

- /Checkpoints/Checkpoint2.ipynb
- /Checkpoints/Checkpoint6.ipynb
parent 8347f674
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Checkpoints/Checkpoint2.ipynb
+ 1
− 1
View file @ 38fe6fe4
%% Cell type:markdown id:75c1fe44 tags:
## Checkpoint 2
### Aim
To write a Python program to read in the coefficients of a quadratic equation and print out the roots for all conditions.
This program demonstrates how to use a simple function, the use of conditional statements, mainly the `if : elif : else:` construct, and how to deal with the various input conditions.
### Task
Write a Python program to read in the three floating point coefficients a, b and c of the quadratic equation:
$$ a x^2 + b x + c = 0$$
then calculate and display the roots using the standard formula for roots of quadratic being,
then calculate and display the roots using the standard formula for roots of a quadratic equation, being :
$$ r_i = \frac{ −b± \sqrt{b^2−4ac}}{2a} $$
The program is expected to deal with conditions of
1. two real roots
2. two complex roots
3. single real root
Your program must contain a function to calculate the discriminant, $d=b^2−4ac$ and must have a main(): method. All variables must be declared inside a function or main().
### Background
In addition to the sections for Checkpoint 1, you will need to have read and studied the following additional sections:
1. [Basic Functions](../CourseNotes/basicfunctions.ipynb)
2. [Conditionals](../CourseNotes/conditionalstatements.ipynb)
Note: the structure of this program is more complex than it initially looks; plan out what you are going to do before you start coding and also check it a bit at a time.
### Checkpoint
- Check your program with the following values:
- a=2.0, b=4.0 and c=−8.0
- a=2.0, b=−4.0 and c=8.0
- a=2.0, b=4.0 and c=2.0
- a=0.0, b=6.0 and c=8.0
Ensure you get sensible answers for the four cases. **Remember when a=0.0 you will get a linear equation with a single root. You should handle this correctly.**
> #### Note:
> You should think about using the `cmath.sqrt` and `math.sqrt` functions to deal explicitly with square roots of complex and real numbers respectively.
- Call a demonstrator, show them your code and run your program with the four sets of inputs above and any additional ones requested by the demonstrator.
- Ensure that the demonstrator finds you in Learn. Please check that your pass/fail (0/1) has been entered correctly (you should be able to see this in Learn).
### Assessment
To pass this checkpoint:
- The program must work correctly for at least three of the above four cases.
- The prompting is sensible.
- The code structure uses a function to calculate the discriminant and has a main(): function.
- The code uses conditional statements `(if : elif : else:)` to test for the various cases.
### Next steps
Once you have completed this checkpoint, please continue with the [Week2](../WeeklyTasks/Week2.ipynb) tasklist.
%% Cell type:code id:883eb6c3 tags:
``` python
```
......
Checkpoints/Checkpoint6.ipynb
+ 1
− 1
View file @ 38fe6fe4
%% Cell type:markdown id:51574bca tags:
## Checkpoint 6
### Aim
To write a Python program to calculate the volume of molecules made up of overlapping spherical atoms; it is also an example in using OOP to simpify the task. Don't get frightened, for this checkpoint you will be given part of the code to help you progress (***this also simulates pretty well how multi-person software development works***).
## Monte Carlo estimate of π
Consider a circle of radius r within a square of side 2r$\times$2r as shown in the figure below.
Consider a circle of radius r within a square of side $2r \times 2r$ as shown in the figure below.
![dynamicpng-width_300.png](dynamicpng-width_300.png)
### Calculation of π by estimation of area.
If we randomly pick N points evenly distributed within the square and $N_{i}$ of them are inside the circle, then an estimate of the area of the circle A is given by the fraction falling inside the circle multiplied by the area of the square, so
$$ A=4r^ 2 \frac{N_i}{N} $$
Noting that A=$\pi r^2$, we get an estimate for the value of π given by,
$$ π≈4 \frac{N_i}{N} $$
So provided the selection of points is evenly distributed over the whole square, than as N becomes large you will get a good estimate for $\pi$. This is the basis of Monte Carlo integration and can be used to estimate the internal area of any closed object.
See an example Python code that uses the random.uniform() function to form random points:
- Code example to calculate $\pi$ [here](../CodeExamples/PiExample.ipynb)
Take a look at this and experiment, you will find that for N>200,000 you get a reasonable estimate for $\pi$, and that as N gets larger the estimate gets better.
This is clearly a "very bad" way to calcuate π, but is very useful for finding the area of overlapping objects in two-dimensions as shown in the figure below; which is otherwise analytically tricky.
![dynamicpng-width_250.png](dynamicpng-width_250.png)
This can also easily be extended to three-dimensions so for example overlapping spheres as shown below.
![dynamicpng-width_2502.png](dynamicpng-width_2502.png)
Here the method is to form a bounding box of dimensions Vb=a×b×c that totally encloses the object. If N points are in the bounding box and N$_i$ are inside the object, then the volume of the object is simply given by
$$ V_0=V_b \frac{N_i}{N} $$
This technique is used in computer simulations of molecules to calculate the volume of molecules that are made up of overlapping spherical atoms.
#### Code Structure
This is rather tricky program but it can be vastly simplied by use of OOP; the suggested structure is
1. Form an Atom class by extending the Vector3d class given by Vector3d Class
2. The Atom class should add an extra internal variable radius and have a boolean method isInside(v) which takes argument v being a Vector3d and returns True if v is inside the atom. (Note the distance(v) method of Vector3d is useful here.)
3. Form a Molecule class that has a list of objects of type Atom and a boolean method isInside(v), which takes argument v being a Vector3d and returns True if v is inside the molecule.
- You also need to consider how to read in the location of the atoms; this is similar to Checkpoint 4.
4. Form a BoundingBox class where the constructor takes a populated Molecule and finds and stores internally the minima / maxima of the box round the molecule (you need 6 floats). This class should have two methods:
- volume() that gets the volume of the box.
- getRandomPoint() which should return a random Vector3d that is randomly distributed inside the box.
(This is the most “messy” bit of the code to write.)
If you have this structure, then what your main program has to do is:
- Ask for the file containing the Atom data and read this into a Molecule.
- Construct a BoundingBox from the Molecule
- Go round a while loop,
- Getting a random Vector3d from the BoundingBox
- Testing if it isInside() the Molecule and incrementing a counter if it is.
- Finally calculating the estimate for the Molecule volume and printing it out.
It is also useful to monitor the progress of the estimate as was done in the example code for estimation of π to check how the integration is going.
#### Format of Data
There are two data files suppied (in the current directory, so they can be read in and executed directly upon here), these being :
- [oxygen.data](oxygen.data) : single oxygen atom
- [water.data](water.data) : water molecule
each consists of a comma delimited CSV file with one atom per line in the format,
Name , r , x , y , z
with comment lines delimited by # being the same format as coded in [Parsing Input](../CourseNotes/ParsingInput.ipynb) and used in [Checkpoint 4](Checkpoint4.ipynb).
#### Computing Task
Write a Python program based on the above structure to:
1. Read in the data file defining the atom locations.
2. Prompt the user for the number of points to run in the calculation.
3. Calculate the volume of the molecule using Monte Carlo integration and output the final value.
4. Optionally plot a graph of the convergence of the algorithm, being the estimate against number of points. The graph should have 100 data points at regular intervals, see Pi example for how to to this easily.
5. Test your code with a single oxygen atom and with water. ***For a pass, only the atom is strictly necessary to run.***
> #### Note
> The volume of a single oxygen atom of radius 1.4 Å is 11.494 Å$^3$. The volume of water is $\approx$ 15.5 Å$^3$. Use these values to check your results before calling a demonstrator.
6. Have another careful read of [final weekly task](../WeeklyTasks/Week6.ipynb).
%% Cell type:code id:9376506d tags:
``` python
""""
A more developed Vector class for simple 3d vectors
"""
import math
import random
class Vector3d(object):
"""
Simple vector 3d class
"""
def __init__(self, x, y, z):
"""
Constructor to form a vector
"""
self.x = float(x)
self.y = float(y)
self.z = float(z)
def __str__(self):
"""
Method to format a vector as a string (implments str())
"""
return "[" + str(self.x) + ", " + str(self.y) + ", " + str(self.z) + "]"
def copy(self):
"""
Method to return a copy of current vector
"""
return Vector3d(self.x,self.y,self.z)
def magnitude(self):
"""
Return the magnitude
"""
return math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
def distance(self,b):
"""
Distance from self to vector b
"""
dx = self.x - b.x
dy = self.y - b.y
dz = self.z - b.z
return math.sqrt(dx*dx + dy*dy + dz*dz)
def normalise(self):
"""
Normalise current vector so its magnitude is 1.0 (unit vector)
"""
mag = self.magnitude()
if mag != 0.0 : # if not already zero
self.x /= mag # dveive throgh by mag
self.y /= mag
self.z /= mag
def add(self,b):
"""
Add a vector b to self and return a new vector.
"""
x = self.x + b.x
y = self.y + b.y
z = self.z + b.z
return Vector3d(x,y,z)
def addTo(self,b):
"""
Add a vector b to self in place.
"""
self.x += b.x
self.y += b.y
self.z += b.z
def dot(self,b):
"""
Form dot product between self and vector b
"""
return self.x*b.x + self.y*b.y + self.z*b.z
def cross(self,b):
"""
Form cross product between self and vector b.
"""
x = self.y*b.z - self.z*b.y
y = self.z*b.x - self.x*b.z
z = self.x*b.y - self.y*b.x
return Vector3d(x,y,z)
def main():
"""
Simple test program to test the vectors
"""
a = Vector3d(1,2,3)
b = Vector3d(4,5,6)
print("Magnitude of a is " + str(a.magnitude()))
c = a.cross(b) # Form cross product
#
# the str() function call the __str__(self) method
print("Cross product of a x b is : " + str(c))
#main()
```
%% Cell type:code id:e0aad9ca tags:
``` python
### Atom class
class Atom(Vector3d):
"""
The Atom class to extends the Vector3d that hold position
"""
def __init__(self,x,y,z,r):
"""
Constructor to form an atom with x,y,z,r parameters
"""
# Your code here
def __str__(self):
"""
Implemnts str()
"""
# Your code here
def isInside(self,v):
"""
isInside method, returns True of v in inside atom
"""
# Your code here
```
%% Cell type:code id:fccc8a4e tags:
``` python
### Molecule class
class Molecule(list):
"""
Molecule, being a list to hold Atoms
"""
def __init__(self,*args):
"""
Constructor to form a Molecule with optional arguments each one assumed to be an Atom
"""
list.__init__(self) # Init underlying list
for a in args: # Append any args to list
self.append(a)
def __str__(self):
""" Format string of all atoms
"""
s = ""
for a in self:
s += str(a) + "\n"
return s
def inInside(self, v):
"""
Method to check if vector is inside moleule.
"""
# Your code here
return False # if here then outside
def fromFile(self,file):
"""
Read Atoms from a file
"""
# Your code here
```
%% Cell type:code id:663fe126 tags:
``` python
### BoundingBox class
class BoundingBox(object):
"""
Class to define a bounding box class round the molecule
"""
def __init__(self, molecule):
"""
The Constructor, which takes a molecule
"""
self.xmin = float("inf") # Initially set limits to +/- infty
self.ymin = float("inf")
self.zmin = float("inf")
self.xmax = float("-inf")
self.ymax = float("-inf")
self.zmax = float("-inf")
for a in molecule: # update limits for each atom in molecule
# Your code here
def volume(self):
"""
Get the volume of the box
"""
# Your code here
def getRandomPoint(self):
"""
Get a random point in the bounding box as a Vector3d.
"""
# Your code here
```
%% Cell type:code id:3e01a08e tags:
``` python
""" Main program to do the simulation
"""
file = open(str(input("File : ")),"r") # Open file
mol = Molecule().fromFile(file) # Read in molecule
print("Molecule is \n" + str(mol)) # Print out info
box = BoundingBox(mol) # Form Bounding box
maxpoint = float(input("Number of points : ")) # Max number of points
plotInterval = maxpoint/100 # Plot internal for monitoring
# Your code here
while p < maxpoint: # Loop counting internal points
# Your code here
...
print("Final estimate is : " + str(estimate))
plt.plot(xData,yData) # Draw graph (with default plot)
plt.ylim(0.9*estimate,1.1*estimate)
plt.title("Estimate of volume for " + str(maxpoint) + " points.")
plt.show()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment