!等方性長方形板のたわみwを計算
! 土木学会『構造力学公式集』をそのまんま書き写す。
! モーメント、せん断力、反力は書き足さないといけない。
! ループm,nは適当に変えるべし。
implicit none
!#--手入力
real(8) :: a=1.72d0    !   m :自由辺になっている辺の長さ(x)
real(8) :: b=3.38d0    !   m :単純支持している辺の長さ(y)
real(8) :: h=240.d-3   !   m :厚さ(z)
real(8) :: P=100.d3    !   N :荷重(z)
real(8) :: u=200.d-3   !   m :載荷幅(x)
real(8) :: v=500.d-3   !   m :載荷幅(y)
real(8) :: xi=1.d0     !   m :載荷中心位置(x)
real(8) :: eta=2.d0    !   m :載荷中心位置(y)
real(8) :: x=1.d0      !   m :求めたいたわみの位置(x)
real(8) :: y=2.d0      !   m :求めたいたわみの位置(y)
real(8) :: E=3.5d9     !  Pa :ヤング率
real(8) :: nu=0.4      !     :ポアソン比
!---
real(8) pi,w0,w1,D,p0,gmm,lmd,alp,beta
real(8) amn,Am,Bm,Cm,Dm,R,F1,F2
integer m,n

pi=2.d0*asin(1.0)!円周率
D=E*h**3/(12*(1-nu**2))!曲げ剛性
p0=P/(u*v)!部分等分布荷重

!---初期化
w0=0.d0

!#--ループm
do m=1,100
   !#--ループn
   do n=1,10
      alp=m*pi/a
      beta=n*pi/b
      lmd=b/a
      !---部分分布荷重
      amn=16.d0*p0/(pi**2.d0*m*n)*sin(m*pi*u/(2.d0*a))*sin(m*pi*xi/a)*sin(n*pi*v/(2.d0*b))*sin(n*pi*eta/b)
      !---Navierの解
      w0=w0+a**4*lmd**4/(D*pi**4)*amn/((lmd**2*m**2+n**2)**2)*sin(alp*x)*sin(beta*y)
   enddo
enddo
print('(A,f17.8)'),'w=',w0
end

