Calculation of the deflection of a stiff beam under gravity

#!/usr/bin/python
##########################################################################
## Program to (try to) calculate the shape of a stiff paper sheet
## deflected due to gravity. Some more info in my blog...

## Coded by Nicolau Leal Werneck <nwerneck@gmail.com> in 2011/04/27

from pylab import *

def pplot(x, args):
    plot(x[:,0], x[:,1], args)

def calc_p(p,dp,phi,ds):
    p[0] = array([0,0])
    for n in range(1,len(p)):
        dp[n-1] = ds*array([cos(phi[n-1]),sin(phi[n-1])])
        p[n] = p[n-1] + dp[n-1]

ion()

Np = 50

## The angles at each point
phi = zeros(Np)

# Set initial state
phi[0]=arctan2(10,0.1)
phi[1:]=0

## Gravitational acceleration
Ag = array([0,-9.8])

## "time interval" between updates
dt = 0.1

## line segment length
ds = 1.0/Np

## Paper linear density
mu = 1 #0.01

## This EI is actually EI / mu
EI = 100

## Upper limit on the number of iterations
Nt = 100000

## A unit vector in teh direction of each segment
dp = zeros((len(phi),2))
## the positions of the weights
p = zeros((len(phi)+1,2))

##
## The main loop

for t in range(Nt):
    ## Calculates the current weights positions.
    calc_p(p,dp,phi,ds)

    ## Variable to test if changes are already too small
    maxchange=0

    ## Iterate over each weight, and calculate its increment
    for s in range(1,len(phi)):
        ## Bending stiffness acceleration
        if s < len(phi)-1:
            Am = EI * ( -p[s-1]+2*p[s]-p[s+1] )  / (2 * ds * mu)
        else:
            Am = EI * ( -p[s-1]+p[s] )  / (2 * ds * mu)

        # Total acceleration
        #At = Am + (Ag if p[s,1]>0 else Ag-1*p[s,1]/mu )
        At = Am + Ag

        dM = dot(At, dot(array([[0,-1],[1,0]]), dp[s]))
        phi[s] += dt*dM
        mm=np.max(np.abs(dM))
        maxchange = mm if mm>maxchange else maxchange

    if maxchange<1e-10:
        break

    if not t%100:
        print t,maxchange

print 'Iterations:', t

calc_p(p,dp,phi,ds)

## Plot the curve
figure(1, figsize=(6,6))
title('Stiff paper sheet bending under gravity')
pplot(p, 'r-o')

plot( [0,10*dp[0,0]], [0,10*dp[0,1]] , 'k--', lw=1)
plot( [0,10*dp[0,0]], [0,10*dp[0,1]] , 'k--', lw=1)
plot([0],[0], 'k-s')

axis('equal')

## The peak of the curve
mm=p[find(p[:,1] == p[:,1].max())[0]]

## Plots an approximation by a parabola
xx = mgrid[0:100.] * p[-1,0]/100;
A = 2 * mm[1]/mm[0]
B = - mm[1]/mm[0]**2
plot(xx, A*xx+B*xx**2, 'g:', alpha=0.5)

savefig('beam3.png', dpi=100)