
To perform a modal analysis we start by constructing the modal force vector from the ground acceleration record.

```python
#Initialise a 2D array to hold effective inertia forces applied to each mass
Peff = np.zeros([N,len(acc)])
#Calculate effective inertia force for each mass
for i in np.arange(N):
    Peff[i,:] = -acc*mArray[i] #(N)

modalP = np.matrix(Phi.transpose()*Peff) # Modal force vector
```

Then we can redefine our piecewise-exact method function.

```python
#Define a function that takes in a time and force vector
def calculateResponse(tVector, fVector):
    '''
    NOTE: All our system parameters are defined globally and not explicitly passed
    into our function as arguments.
    '''

    nPoints = len(tVector) #Number of data points
    delta_t = tVector[1]-tVector[0]

    #Constants (ref. Dynamics of Structures by Humar)
    A = math.e**(-xi*wn*delta_t)*((xi/(math.sqrt(1-(xi**2))))*
        np.sin(wd*delta_t)+np.cos(wd*delta_t))

    B = math.e**(-xi*wn*delta_t)*((1/wd*np.sin(wd*delta_t)))

    C = (1/k)*(((2*xi)/(wn*delta_t)) + math.e**(-xi*wn*delta_t)*
        ((((1-2*(xi**2))/(wd*delta_t))-(xi/(math.sqrt(1-(xi**2)))))*
        math.sin(wd*delta_t)-(1+((2*xi)/(wn*delta_t)))*np.cos(wd*delta_t)))

    D = (1/k)*(1-((2*xi)/(wn*delta_t))+math.e**(-xi*wn*delta_t)*
        (((2*(xi**2)-1)/(wd*delta_t))*np.sin(wd*delta_t)+
        ((2*xi)/(wn*delta_t))*np.cos(wd*delta_t)))

    A1 = -math.e**(-xi*wn*delta_t)*((wn/(math.sqrt(1-(xi**2))))*
         np.sin(wd*delta_t))

    B1 = math.e**(-xi*wn*delta_t)*(np.cos(wd*delta_t)-((xi)/
         (math.sqrt(1-(xi**2))))*np.sin(wd*delta_t))

    C1 =(1/k)*(-(1/delta_t)+math.e**(-xi*wn*delta_t)*
        ((((wn)/(math.sqrt(1-(xi**2))))+((xi)/(delta_t*
        math.sqrt(1-(xi**2)))))*np.sin(wd*delta_t)+(1/delta_t)*np.cos(wd*delta_t)))

    D1 =(1/k)*((1/delta_t)-(math.e**(-xi*wn*delta_t)/delta_t)*
               ((xi/(math.sqrt(1-(xi**2))))*np.sin(wd*delta_t)+np.cos(wd*delta_t)))

    #INITIAL CONDITIONS
    x0=0
    v0=0

    pos = [] #Initialise an array to hold calculated position
    vel = [] #Initialise an array to hold calculated velocity

    #MAIN LOOP (to 'step through' solution)
    for n, timestep in enumerate(tVector):
        #note that enumerate adds a counter to an iterable (n)
        timestep = round(timestep,3) #The current time

        #Update the force values
        if n<nPoints-1:
            Fn = fVector[n]
            Fnp1 = fVector[n+1]
        else:
            #Pickup when we're at the end of the force vector (no more +1 values)
            Fn = fVector[n]
            Fnp1 = 0

        #Calculate the position and velocity at the end of the timestep
        curPos = (A*x0) + (B*v0) + (C*Fn) + (D*Fnp1)
        curVel = (A1*x0) + (B1*v0)+ (C1*Fn) + (D1*Fnp1)

        #Append (save) the position and velocity into the vectors defined above
        pos.append(curPos)
        vel.append(curVel)

        #Update initial conditions for the next loop iteration
        x0 = curPos;
        v0 = curVel;

    return pos, vel #Return the calculated position and velocity out of the function
```

### Calculate modal and structural responses

We use a `for` loop to call our `calculateResponse()` function on each mode and save the modal response. Once all of the modal responses have been calculated we can recombine them in a single line using the modal matrix.

```python
modalResponses = np.zeros([N,len(acc)]) #Initialise a 2D array to hold modal responses
for i in np.arange(N):
    #Redefine modal parameters used in the calculateResponse function

    wn = Omega[i]
    k = wn**2
    xi = xiArray[i]/100
    wd = wn*math.sqrt(1-xi**2)

    modalLoad = np.ravel(modalP[i,:]) #Flatten modal load (matrix -> array)
    pos, vel = calculateResponse(time, modalLoad) #Call the function
    modalResponses[i,:] = pos #Save the response

#Recombine modal responses to obtain structural response
U = Phi*modalResponses
```

Next we plot the modal responses. As we might expect, the structural response is governed by the first mode. This first mode sway response typically dominates for shear buildings experiencing ground motion.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,2])

for i in np.arange(N):
    axes.plot(time, modalResponses[i,:], label=f"Mode {i+1}, {round(Omega[i],2)} rads/sec")

axes.set_xlim([0,max(time)])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Modal amplitude')
axes.set_title('Modal responses')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/ground-motion-response-using-modal-superposition/img1.jpg"
	alt="Modal responses. | EngineeringSkills.com"
	width="1000"
	height="688"
	caption="Fig 1. Modal responses."
/>

Plotting the structural response confirms the dominance of the first mode.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,2])

for i in np.arange(N):
    axes.plot(time, np.ravel(U[i,:]), label=f"Mass/Floor {i+1}")

axes.set_xlim([0,max(time)])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Structural responses')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/ground-motion-response-using-modal-superposition/img2.jpg"
	alt="Structural response of each mass/floor. | EngineeringSkills.com"
	width="1000"
	height="690"
	caption="Fig 2. Structural response of each mass/floor."
/>


Please log in or enroll to continue

Please log in or enroll to continue