fcoord = [0.0 for i in range(10)]
fmag = [0.0 for i in range(10)]
mcoord = [0.0 for i in range(10)]
mmag = [0.0 for i in range(10)]
dcoord1 = [0.0 for i in range(10)]
dcoord2 = [0.0 for i in range(10)]
dint = [0.0 for i in range(10)]
slf = [0.0 for i in range(10)]
sld = [0.0 for i in range(10)]
slm = [0.0 for i in range(10)]
df = [0.0 for i in range(10)]
dd = [0.0 for i in range(10)]
dm = [0.0 for i in range(10)]
sl = [0.0 for i in range(100)]
d = [0.0 for i in range(100)]
ddist = [0.0 for i in range(10)]
load = [0.0 for i in range(10)]
slfx = 0.0
sldx = 0.0 
slmx = 0.0
dfx = 0.0
ddx = 0.0
dmx = 0.0


if 0:
    # youngs modulus
    e = 200.0e9
    iner = (0.06)*(0.28)**3/12
    len = 12.
    nf = 1
    nd = 0
    nm = 1
    num = 25
    rcoord1 = 0
    rcoord2 = 8.
    fcoord[0] = 12.
    fmag[0] = 8500.
    mcoord[0] = 4.
    mmag[0] = 8000.
    #dcoord1[0] = 0
    #dcoord2[0] = 0
    #dmag[0] = 0

    [12,8500]
    [4,800]
    r1 = float(5250)
    r2 = float(-13750)
if 1:
    e = 30e6
    iner = 3.375
    len = 100
    nf = 1
    nd = 1
    nm = 1
    num = 50
    rcoord1 = 20.
    rcoord2 = 50.
    fcoord[0] = 0.
    fmag[0] = 3000.
    mcoord[0] = 100.
    mmag[0] = -8000.
    dcoord1[0] = 20.
    dcoord2[0] = 50.
    dint[0] = 125.

def isec(b, s, t, d):
    h = d - 2*s
    y = d/2
    A = b*d - h*(b-t)
    I = (b*d**3 - h**3*(b-t))/12
    print I,A
    return I

def tsec(b, s, t, d):
    h = d - s
    y = d - (d**2*t + s*s*(b-t))/(2*(b*s+h*t))
    I = (t*y**3 + b*(d-y)**3 - (b-d)*(d-y-s)**3)/3
    print y, I
    return I/1e12

isec(200,9,9,200)

if 1: # kitchen beam
    e = 200e9
    iner = (1.2**4)*2140/1e8#tsec(200, 6, 6, 200)
    len = 3.
    nf = 3
    nd = 1
    nm = 0
    num = 50
    rcoord1 = 0.0
    rcoord2 = 3.0
    fcoord[0] = 0.0
    fmag[0] = 000.
    fcoord[1] = 1.45
    fmag[1] = 17000.
    fcoord[2] = 3.0
    fmag[2] = 000.
    #mcoord[0] = 100.
    #mmag[0] = -8000.
    dcoord1[0] = 0.
    dcoord2[0] = 3.3
    dint[0] = 17000

print 'the magnitude and direction of the pin reaction forces'

fsum = 0.0
msum = 0.0
for i in range(nd):
    ddist[i] = dcoord2[i] - dcoord1[i]
    load[i] = dint[i]*ddist[i]
    fsum = load[i] + fsum
    msum = (((0.5*ddist[i] + dcoord1[i]) - rcoord1) * load[i]) + msum

for i in range(nf):
    fsum += fmag[i]
    msum += ((fcoord[i] - rcoord1)*fmag[i])

for i in range(nm):
    msum += mmag[i]

r2 = -(msum/(rcoord2-rcoord1))
r1 = -(fsum + r2)

print 'r1 = ', r1, 'r2 = ', r2

print 'calculating the largest number of either force, distributed loads, or moments'

big = max(nf, nd, nm)

dx = float(len)/num
print 'len = ', len, 'dx = ' , dx
j = 0
for x in [rcoord1, rcoord2] + [i*dx for i in range(num+1)]:
    for i in range(big):
        vv1 = max(0.,x - fcoord[i])
        vv2 = max(0.,x - dcoord1[i])
        vv3 = max(0.,x - dcoord2[i])
        vv4 = max(0.,x - rcoord1)
        vv5 = max(0.,x - rcoord2)
        vv6 = max(0.,x - mcoord[i])
        
        slf[i] = fmag[i]*(vv1**2)/2
        df[i] = fmag[i]*(vv1**3)/6
        sld[i] = (dint[i]*(vv2**3) - dint[i]*(vv3**3))/6
        dd[i] = (dint[i]*(vv2**4) - dint[i]*(vv3**4))/24
        slm[i] = mmag[i]*vv6
        dm[i] = mmag[i]*(vv6**2)/2
        
    slr1 = r1 * (vv4**2)/2
    slr2 = r2 * (vv5**2)/2
    dr1 = r1 * (vv4**3)/6
    dr2 = r2 * (vv5**3)/6

    for i in range(big):
        slfx += slf[i]
        sldx += sld[i]
        slmx += slm[i]
        dfx += df[i]
        ddx += dd[i]
        dmx += dm[i]
    
        sl[j] = slfx + sldx + slmx + slr1 + slr2
        d[j] = dfx + ddx + dmx + dr1 + dr2
    
    slfx = 0.0
    sldx = 0.0
    slmx = 0.0
    dfx = 0.0
    ddx = 0.0
    dmx = 0.0
    j += 1

c1 = (d[1] - d[0])/(rcoord1 - rcoord2)
c2 = (-d[0] - (c1*rcoord1))

x = 0.0

for i in range(2, j):
    sl[i] = (sl[i] + c1)/(e*iner)
    d[i] = (d[i] + (c1*x) + c2) / (e*iner)
    #print sl[i], d[i]
    x += dx

x = 0.0
for i in range(2, j):
    print i-2, x, "%.3g %.3gmm" % (sl[i], d[i]*1e3)
    x += dx