Solve Linear Equations in Three Variables - Using Python, Numpy and Sympy , With the Given Code

 Hi,

In the series of posts about Python for Civil Engineers, I have come with something from linear algebra.

While dealing with statics problems, such as finding unknown forces in space truss members, one has to deal with the system of linear equations in three variables. 

One can use the 'Sympy' and 'Numpy' libraries in Python can get a beautiful solution code for solving a system of linear equations in three variables. Computer programming is very much easy in Python, and you definitely need to try it. The code given below is separated into two parts: 1st part is the solution of the same problem using Sympy, and 2nd part gives solution code using Numpy library. 

The code is written specifically to solve the following three linear equations, but you can change the numbers accordingly.

x + 2y + 3z = 8

4x + 2y + 4z = 5

4x + 2y + 5z = 1

Simply change the coefficient for any other system. 

Code:

import numpy as np
import sympy as sp
from sympy import Eq, solve_linear_system, Matrix

#  Using sympy

x,y,z = sp.symbols('x y z')
row1 = [1, 2, 3, 8]
row2 = [4, 2, 4, 5]
row3 = [4, 2, 5, 1]
system = Matrix((row1, row2, row3))
print(system)
print(solve_linear_system(system,x,y,z))

print('___________________________________ ')

# Using numpy
row1 = np.array([1,2,3])
row2 = np.array([4,2,4])
row3 = np.array([4,2,5])
constants = np.array([8,5,1])
system = np.array((row1, row2, row3))
print(system)
xyz= np.linalg.solve(system, constants)
print(xyz)

Thanks!

Python Code to Find out the Dot Product of Vectors and Angle between them (Using end Point Co-ordinates too)

Hi, 

Sometimes you are given with the endpoints' co-ordinates of vectors, and sometimes you are directly given the vectors. For both conditions, the code given below will work. You have to input them as and when asked. End results will give the vectors, their magnitudes, their dot product and the angle between the two vectors. 

Python is a computer programming language that can be used to develip softwares, apps, and websites etc. The following python code provides a short solution to the Vector Algebra. 

Code: 

import numpy as np

import math

 a1 = input('For vector1, co-ordinates given? (y/n): ')

a2 = input('For vector2, co-ordinates given? (y/n): ')
u = input('unit: ')
if a1 == 'y':
    x11= float(input('enter co-ordinates of starting point, x11=  '))
    y11 = float(input('y11= '))
    z11 = float(input('z11 = '))
    V11 = np.array([x11, y11, z11])

    x21 = float( input('enter co-ordinates of end point, x21=  '))
    y21 = float(input('y21= '))
    z21 = float(input('z21 = '))
    V21 = np.array([x21, y21, z21])

    print(f'Starting point, V11 = {V11}')
    print(f'Ending point, V21 = {V21}')
    v1 = V21 - V11
    print(f'Vector1, V1 = V21 - V11 = {v1}')
    v1m = np.sqrt((x21-x11) ** 2 + (y21-y11) ** 2 + (z21-z11) ** 2)
    print(f" Magnitude of v1, v1m = {v1m} {u}")
else:
   print('enter the elements in vector1 ')
   x1 = float(input())
   y1 = float(input())
   z1 = float(input())
   v1 = [x1,y1,z1]
   print(f'v1 =({x1}i + {y1}j + {z1}k) {u}' '\n')
   v1m = np.sqrt(x1**2 + y1**2 + z1**2)
   print(f" Magnitude of v1, v1m = √(({x1})^2 + ({y1})^2 + ({z1})^2 = {v1m} {u}" '\n')

if a2 == 'y':
    x12= float( input('enter co-ordinates of starting point, x21 =  '))
    y12 = float(input('y12= '))
    z12 = float(input('z12 = '))
    V12 = np.array([x12, y12, z12])

    x22 = float( input('enter co-ordinates of end point, x22 =  '))
    y22 = float(input('y22= '))
    z22 = float(input('z22 = '))
    V22 = np.array([x22, y22, z22])
    print(f'Starting point, V12 = {V12}')
    print(f'Ending point, V22 = {V22}')
    v2 = V22 - V12
    print(f'Vector2, V2 = V22 - V12 = {v2}')
    v2m = np.sqrt((x22 - x12) ** 2 + (y22 - y12) ** 2 + (z22 - z12) ** 2)
    print(f" Magnitude of v2, v2m = {v2m} {u}")

else:
    print('enter the elements in the 2nd vector: ')
    x2 = float(input())
    y2 = float(input())
    z2  = float(input())
    v2 = [x2,y2,z2]
    print(f'v2 = ({x2}i + {y2}j + {z2}k) {u}' '\n')
    v2m = np.sqrt(x2**2 + y2**2 + z2**2)
    print(f" Magnitude of v2, v2m = √(({x2})^2 + ({y2})^2 + ({z2})^2 = {v2m} {u}" '\n')

d = float(np.dot(v1,v2))
e = d/(v1m*v2m)
a = math.acos(e)
k = math.degrees(a)

print(f'''Taking Dot product, v1.v2 = v1m * v2m * cosθ 

      => cosθ = v1.v2/(v1m*v2m) = {d}/({v1m}*{v2m}) = {e} 

      hence, Angle between the two vectors, θ = acos({e}) = {k}
''')

Python Code - Discharge in Rectangular Channel using Manning's Formula

 Hi,

Next in the series 'Python for Civil Engineers' is the code for Manning's formula for the rectangular channel.

One has to provide the needed inputs to get the discharge, in the given rectangular channel, as the output. The Manning's formula for discharge calculations can be used for various channels such as rectangular, trapezoidal etc. The Python code in this post is written for the rectangular channel, in particular. Inputs needed are Manning's n, slope s, bottom width b, and water depth h. Also, mention the units when asked. 

def discharge():
    n = float(input("n = Manning's value = "))
    S = float(input("S = Channel's slope = "))
    u = input("unit of width/height (m or ft)= ")
    b = float(input("b = bottom width = "))
    h = float(input("h = water depth = "))
    k = 1
    u1 = "Cumecs"
    if u.lower() == "ft":
        u1 = "cfs"
    if u.lower() == "ft":
        k = 1.49

    Q = (k/n) * (b * h) * (b * h / (2.0 * h + b)) ** (2.0 / 3.0) * S ** (0.5)
    Q = round(Q,3)
    print(f' Q = {Q} {u1}')
    return Q

print('''
    Manning's Formula for discharge in rectangular channel Q = (k * A * R^(2/3) * S^0.5) / n
        Q = Flow Rate, (in cfs, or m)
        k = 1.49 for discharge in cfs, and 1 for in cumecs
        A = Area of section = b*h , b = bottom width, h = water depth
        R = hydraulic Radius = Area/Wetted Perimeter = bh/(2h+b)
        n = Manning's Roughness coefficient
        ''')

print(f'Inputs needed to calculate the discharge Q')
discharge()

Thanks!

Vertical Equal Tangent Parabolic Curve - Solution Code in Python (Computer Program)

 Hi,

Now solving the vertical parabolic curve for the stations and elevations at the full stations is easier with the coding. A computer program that I have coded for you runs on Pycharm (Python). Pycharm is a software that is used for offline coding on your personal computer. If you have a good Internet connection, you can produce the code on online platforms such as Jupyter. 

 You have to put the value of the grades, length of the curve, station PVC and Elevation of PVC. It will calculate all other stations and elevations for you. 

Code: 

import math

g1 = input('Grade G1 (%) = ')
a = float(g1)/100
g2 = input('Grade G2 (%) = ')
b = float(g2)/100
Length = input('Length of the vertical curve = ')
L = float(Length)
u = input('unit = ')
A = float((b - a)/L)
print('Rate of change of grade = k = (G2-G1)/L =', A)
print('Let x be the distance from PVC to a point on the curve,'
      ' and Y be the respective elevation.')
c = float(input('Station(PVC) = Xo = ',))
d = float(input('Elevation at PVC = Yo = ',))
print(f"Equation of a parabolic equal tangent curve, "
      f"Y = Yo + G1x + kx^2/2 = {d} + {a}x + {A}x^2/2     ..eq1")
x = 100*((1+c//100) - c/100)
X = c + x,
print('1st Full Station is', X, 'at x =', x, u, 'from PVC')
Y = d + a*x + A*(x**2)/2
print('Putting the value of x in eq1, Elevation, Y= ', Y, u)
x = x + 100
while x < L:
    X = c + x
    Y = d + a * x + A * (x ** 2) / 2
    print('Similarly, next Full Station is', X, 'at x =', x, u, 'from PVC')
    print('Elevation, Y = ', Y, u)
    x = x + 100
X = c + L
print('End Station, EVC = PVC + L =', X, u, 'at x = L = ', L, u)
x = L
Y = d + a*x + A*(x**2)/2,
print('Put x in eq1, Elev EVC = ', Y, u)

Output Sample: 

Grade G1 (%) = 3
Grade G2 (%) = -2
Length of the vertical curve = 1000
unit = m
Rate of change of grade = k = (G2-G1)/L = -5e-05
Let x be the distance from PVC to a point on the curve, and Y be the respective elevation.
Station(PVC) = Xo = 340
Elevation at PVC = Yo = 100
Equation of a parabolic equal tangent curve, Y = Yo + G1x + kx^2/2 = 100.0 + 0.03x + -5e-05x^2/2     ..eq1
1st Full Station is (400.0,) at x = 60.00000000000001 m from PVC
Putting the value of x in eq1, Elevation, Y=  101.71 m
Similarly, next Full Station is 500.0 at x = 160.0 m from PVC
Elevation, Y =  104.16 m
Similarly, next Full Station is 600.0 at x = 260.0 m from PVC
Elevation, Y =  106.11 m
Similarly, next Full Station is 700.0 at x = 360.0 m from PVC
Elevation, Y =  107.56 m
Similarly, next Full Station is 800.0 at x = 460.0 m from PVC
Elevation, Y =  108.50999999999999 m
Similarly, next Full Station is 900.0 at x = 560.0 m from PVC
Elevation, Y =  108.96 m
Similarly, next Full Station is 1000.0 at x = 660.0 m from PVC
Elevation, Y =  108.91 m
Similarly, next Full Station is 1100.0 at x = 760.0 m from PVC
Elevation, Y =  108.36 m
Similarly, next Full Station is 1200.0 at x = 860.0 m from PVC
Elevation, Y =  107.31 m
Similarly, next Full Station is 1300.0 at x = 960.0 m from PVC
Elevation, Y =  105.76 m
End Station, EVC = PVC + L = 1340.0 m at x = L =  1000.0 m
Put x in eq1, Elev EVC =  (105.0,) m

Thanks!

Check Points of Buying a Flat/House

Hi,
Building a house is a mutual work of the builder, government and the buyers. A civil engineer plays its role in designing and constructing the house, and once he does his job, it is the buyer who has the last say if he wants to buy the house or not. If the house looks good, a buyer would love to buy the house, but there are other certain things that must be checked before buying the house from a builder. This blog post is about those important things that I have learned after buying my own flat but should have known them before.

 Check if the builder has got the project registered in the RERA website of the concerned state. For example, I bought a flat in Punjab, I checked the Punjab RERA website to see the project name under the registered projects list. It is very important that you do it because this will probably save you from any kind of possible cheating and fraud.

Things to check in the builder's RERA profile:

(A)
See if the project registration period has not expired. If it is expired, then you must avoid buying the flat.

(B)
See if there is a specified format of 'Agreement of Sale' in the uploaded documents section. You must insist to use this agreement format to finalise a deal with the builder. It has been seen that the builder asks you to sign the agreement as per his own terms and clauses which most of the times would be partial towards you. On the other hand, the format given by RERA will have the clauses which will give you all the rights that RERA has to offer to the buyers.
If you sign an agreement different from the one given by the RERA it is highly likely that builder would harass you in some way, like forfeiting your booking amount in case of cancellation of the deal, asking you for unfair maintenance charges, having partial interest rates for compensation in case of delay in possession of the flat etc.

(C)
If the project is shown as 'completed' in the RERA website then check for the completion and occupancy certificates. They would be uploaded in the website and hence you must check it. Remember that even if you get the registry papers of the flat you are still at risk. Until and unless you have the occupancy and completion certificates you are at the risk of losing your flat. You want to know why!? Please read about Campa Cola building in Mumbai.

(D)

 Don't give more than 10% of the money as the booking amount, this is the amount specified by RERA. If the builder asks you for more, please ask him to read the RERA rules.

(E)

 Ask the builder for the layout plan,  and as-built-plan and layout of the society.

(F)

 Check if the builder has all the clearance certificates and NOCs uploaded on the RERA account. You can see all these papers for free by signing up on the RERA website.

Once everything is okay, only then you must invest your money.

Thanks! 

All the basic Soil Water and Air Relations and Formulae, Three phase diagram II Soil Mechanics

Basic Soil and Water Relationships.
 Soil is considered a three-phase material, made of soil, water and air. The basic formulae used are written below. 

Sl No	Relationship in mass density	Relationship in unit weight
1	n = e/(1+e)	n = e/(1+e)
2	e = n/(1-n)	e = n/(1-n)
3	na = n ac	na = n ac
4	P= (G+Se)Pw/(1+e)	Y= (G+Se) Y w/(1+e)
5	Pd = GPw/(1+e)	Yd = G Yw/(1+e)
6	Psat = (G+e)Pw/(1+e)	Ysat = (G+e) Yw/(1+e)
7	P’ = (G-1)Pw/(1+e)	Y’ = (G-1) Y w/(1+e)
8	e = wG/s	e = wG/s
9	Pd = P/(1+w)	Yd = Y/(1+w)
10	Pd = (1-na)G Pw/(1+wG)	Yd = (1-na)G Y w/(1+wG)

Notations used: 
n = porosity, 
e = void ratio, 
na = air content,
P = Density,  Pd= dry density, 
Psat = saturated density, P' = submerged density
Y = Unit weight density, Yd = dry unit weight, 
Ysat = Saturated Unit weight, Y' = submerged unit weight
G = Specific Gravity
S = Saturation
w = water content

Thanks!

How To Operate a Total Station - Steps - All components explained

Hi,

Total Station(TS) is a surveying instrument used to find the direction angles and distances of various lines. In this post, you will learn how to operate a Total Station. But first, let's know the components of a TS. 

Components.
These images show the components of TS. The carrying handle can be fixed using to handle tightening screws. The sighting collimator is used for the rough sight of the prism. The Eyepiece is used for sighting the prism. Together with the eyepiece is telescope holding grip which is used to transit the telescope. 



All these are shown in the image below. Telescope holding grip is used for transiting the telescope, and telescope focusing knob is used to clear the image on the diaphragm.



Then there is a vertical tangent screw used for the small vertical movement of the TS so as to get the clear and straight line of sight so that the crosshairs intersect the prism. Also, there is a vertical motion clamp used to fix/lock the position of the vertical tangent screw.



Similarly, there is a horizontal tangent screw and a horizontal motion clamp, which is used to control the horizontal rotation of the TS. 
There are two display units, and the keys adjacent to them. The optical plummet screw is used to sight the ground and is used to level the TS.  



There are three levelling screws in the tri-bach which is fixed using the tri-bach fixing lever. 



When you power on the Display unit, it gets into its default mode - i.e.angle measurement mode. Display unit uses dot matrix LCD, which has four lines(rows) and 20 characters per line. The upper three lines display the measured data, and bottom line displays the softkeys' functions which change with the measuring mode.


The image shown above shows all the keys with their respective functions. 

Now before we take any measurements, it is important to do the centering of the TS and the prism. It is done as explained below.

Step.1 Level the Total Station(TS)

Put the Total Station on the Tripod Stand. The approximately level the lower bubble by arranging the three legs of the tripod stand.


Then note that there is a horizontal bubble tube, in which the bubble will not be at centre if the total station is not levelled. 


Bring this bubble tube parallel with the imaginary line joining two of the three screws. And then rotate both of them together either inwards or outwards. Note the movement of the bubble and rotate the screws until the bubble comes to the centre. 

Now rotate the total station such that the bubble tube is perpendicular to that imaginary line between the two screws. When you do so, the bubble may move from the central position, so now rotate the third screw until it again comes to the centre. Now bring the tube to the initial position and again check the bubble and repeat the process until the bubble is at the centre in both of the positions. 
Now your Total Station is Levelled. 

Step.2 Put the Stand and Prism at the station of interest and Level it

There is a stand and a prism which is sighted from the TS. Fix the prism on the stand, and note that there a bubble on the stand too, so you have to level the stand too. 


Simply fix the vertical position of the stand until the bubble is centred.  Now sight the TS through the sighting collimator. 




Now that centring is done, let's take the measurements to find the horizontal distance between two points.

Horizontal Distances, slope distance, co-ordinates and angle Measurement using Total Station. 

So here we have two points, one - where the instrument is set up and the other where the prism is. Use the sighting collimator to do the approximate bisection of the prism.
 Then sight through the eyepiece and use the focussing screw to bisect the prism with the crosshairs. Then tighten up the vertical tangent screw using the vertical tangent clamp. 
Use the eyepiece to get a clear prism image. Once it is done, click on the distance measurement mode key.

Display Unit of Total station showing the Slope Distance(SD) 

Similarly, when you click the co-ordinate key on, it will show the coordinates (x,y,z) of the station. And when it is in angle measurement mode, it will display horizontal angle with the right face(HR) and vertical angle(V.) 

Thanks!