###########################################################################
# program: plate_corners.py
# author: Tom Irvine
# version: 1.0
# date: July 21, 2011
# description:  This program calculates the fundamental bending frequency
#               of a plate supported at each corner.
###########################################################################

from tompy import material
from tompy import enter_float
from tompy import GetInteger2
from tompy import MatrixMax

from math import pi

from numpy import sqrt,zeros,sin,cos,arange,meshgrid

from numpy.random import rand

import matplotlib.pyplot as plt


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FixedLocator, FormatStrFormatter


###########################################################################

E,rho,mu=material()

fmin=1.0e+90  
#
print(' ')  
print(' Enter thickness (in)')  
h=enter_float()  
# 
print(' ')  
print(' Enter length (in)')  
LLL=enter_float()    
# 
print(' ')  
print(' Enter width (in)')  
WWW=enter_float()
#
a=max(LLL,WWW)
b=min(LLL,WWW)    
# 
added_mass=0 
print(' ')  
print(' Add nonstructural mass?  1=yes  2=no')  
am=GetInteger2()

if(am==1): 
   print(' ')  
   print(' Enter mass (lbm) ')  
   added_mass=enter_float()  
   added_mass=added_mass/386  
  
# 
volume = a*b*h  
rho=rho + (added_mass/volume) 
# 
total_mass=rho*volume*386 
print '\n total mass = %8.4g lbm \n' %total_mass 
 
D=E*(h**3)/(12*(1-mu**2))  
DR=D/(rho*h/2)  

print ' Plate Stiffness Factor = %8.4g lbf in\n' %D  
# 
iflag=0

num=40  
nnn=num 
nnnp1=nnn+1

dx=a/num  
dy=b/num  
pia=pi/a  
pib=pi/b  
#
ad2=a/2
bd2=b/2 
#
knum=1000  
alpha_r=0  
beta_r=0 
gamma_r=0 
#
fn1=1.13*sqrt(D/(rho*h))/a**2 
fn2=1.13*sqrt(D/(rho*h))/b**2
  
flow=min(fn1,fn2) 
fup=max(fn1,fn2) 
#
print(' ') 
print ' Expected:  %8.4g < fn < %8.4g Hz ' %(flow,fup) 

#
print(' ') 
print(' Calculating natural frequency.... ') 
#
c1=2*mu 
c2=2*(1-mu) 
#
x=zeros(nnnp1,'f')
y=zeros(nnnp1,'f')

for i in range (0,nnnp1):
    x[i]=(i*dx-(ad2)) 
    y[i]=(i*dy-(bd2))
    
zmax=3 
#
print(' ') 
for k in range (1,knum):
    alpha= 0.8*rand()+0.2 
    beta = 0.8*rand()+0.2 
#

    if(iflag==1 and rand() > 0.5 ):
        an=(float(k)/float(knum))**(0.12)
        dt=2*(1-an)
        alpha=alpha_r*(an+dt*rand()) 
        beta =beta_r*(an+dt*rand())        

    gamma= min(alpha,beta)*rand() 
#    
    if(abs(a-b)<a/100):
        alpha=beta 
    
#
    tsum=alpha+beta+gamma
#    
    alpha/=tsum 
    beta/=tsum 
    gamma/=tsum 
#
    pia_alpha=pia*alpha 
    pib_beta =pib*beta 
    pia2_alpha=(pia**2)*alpha 
    pib2_beta =(pib**2)*beta   
    pia_alpha_pib_beta=pia_alpha*pib_beta
#
    T=0 
    V=0


    x_theta=pia*x
    y_theta =pib*y

    cx=cos(x_theta)
    sx=sin(y_theta)
    
    cy=cos(x_theta)
    sy=sin(y_theta)      
      

    gamma_sx=gamma*sx 
    one_plus_gamma_sx=1+gamma_sx 

    gamma_sy=gamma*sy 
    one_plus_gamma_sy=1+gamma_sy 
      
    for i in range (0,nnnp1):
          
        for j in range (0,nnnp1):
            
            Z=alpha*cx[i] + beta*cy[j] + gamma*cx[i]*cy[j] 
#
            dZdx  =  -pia_alpha*(one_plus_gamma_sy[j])*sx[i] 
            d2Zdx2=  -pia2_alpha*(one_plus_gamma_sy[j])*cx[i]    
#
            dZdy  =  -pib_beta*(one_plus_gamma_sx[i])*sy[j] 
            d2Zdy2=  -pib2_beta*(one_plus_gamma_sx[i])*cy[j]   
#    
            d2Zdxdy=pia_alpha_pib_beta*gamma*sy[j]*sx[i] 
#
            V+=(d2Zdx2)**2 + (d2Zdy2)**2 + c1*(d2Zdx2)*(d2Zdy2) + c2*(d2Zdxdy)**2 
#            
            T+= Z**2 
  
    
    V*=D/2 
    T*=rho*h/2 
    omega=sqrt(V/T) 
    fn=omega/(2*pi) 
    if(fn<fmin and fn>= flow):
        iflag=1
        fmin=fn 
        fnr=fn
        alpha_r=alpha 
        beta_r =beta 
        gamma_r=gamma 
        print '%d %7.4g Hz  alpha=%6.3g  beta=%6.3g  gamma=%6.3g ' %(k,fn,alpha_r,beta_r,gamma_r) 
    

xx=zeros(41,'f') 
yy=zeros(41,'f') 
zzr=zeros((41,41),'f')
#
for i in range (0,41):
    xx[i]=(i*dx-(a/2)) 
    cx=cos(xx[i]*pia) 
    sx=sin(xx[i]*pia) 
    for j in range (0,41):
        yy[j]=(j*dy-(b/2)) 
        sy=sin(yy[j]) 
        cy=cos(yy[j]*pib) 
        zzr[i][j]=alpha_r*cx + beta_r*cy + gamma_r*cx*cy 
               

################################################################################

print " "
print " view plot"

zmax=MatrixMax(zzr)
             
aaa=max(max(xx),max(yy))
xmin=-aaa 
ymin=-aaa 
zmin=0 
xmax=aaa 
ymax=aaa 

dx=xmax/2
dy=ymax/2


fig = plt.figure(1)
ax = fig.gca(projection='3d')

X = arange(xmin, xmax, dx)
Y = arange(ymin, ymax, dy)

X, Y = meshgrid(xx, yy)
Z = zzr

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)


ax.set_xlim3d([-aaa, aaa])
ax.set_ylim3d([-aaa, aaa])


#  ax.w_zaxis.set_major_locator(LinearLocator(10))
#  ax.w_zaxis.set_major_formatter(FormatStrFormatter('%.03f'))

fig.colorbar(surf, shrink=0.5, aspect=5)

string = "fn= %7.3g Hz" % fnr

plt.title(string)

plt.show()