!断面の軸方向応力分布の理論値sigma=M/I・yとsigma=0.のデータをgnuplot用に出力する。

implicit none
!#--入力
real(8) :: p=100.d0                ! N    :荷重
real(8) :: ell=500.d-3             ! m    :スパン
real(8) :: h=10.d-3                ! m    :桁高
integer :: nh=10                   !　    :高さ方向分割数
real(8) :: b=10.d-3                ! m    :幅
integer :: nb=10                   !　    :幅方向分割数
real(8) :: length=0.5d-3           ! m    :1要素の長さ
!---
real(8)    mtan                    ! N・m :固定端の曲げモーメント力
real(8)    mhobotan                ! N・m :固定端の隣の曲げモーメント力
real(8)    mnaka                   ! N・m :梁中央の曲げモーメント力
real(8)    moi                     ! m^4  :断面2次モーメント
real(8)    pi                      ! rad  :円周率
real(8),allocatable :: szz(:,:,:)  ! Pa   :軸方向応力
integer i,j
allocate(szz(0:nb,0:nh,3))
pi=2.d0*asin(1.d0)
!#--出力ファイル名
open(11,file='zero.d')
open(12,file='kotei.d')
open(13,file='hobokotei.d')
open(14,file='naka.d')
!---

!#--断面2次モーメント(電卓等で計算)
!moi=
!moi=(b*h**3.d0)/12.d0 !矩形のみ
moi=pi*(b**2.d0*h**2.d0)/64.d0 !円形のみ

!---応力計算
mtan=-p*ell
mhobotan=-p*(ell-length)
mnaka=-p*ell/2.
do j=0,nh
   do i=0,nb
      szz(i,j,1)=mtan/moi*(h*j/nh-h/2.d0)/(1.d6) !MPaに直した
      szz(i,j,2)=mhobotan/moi*(h*j/nh-h/2.d0)/(1.d6) !MPaに直した
      szz(i,j,3)=mnaka/moi*(h*j/nh-h/2.d0)/(1.d6) !MPaに直した
   enddo
end do
        
!---出力(長さの単位はmmに直した)
do j=0,nh
   do i=0,nb
      write(11,*) (b*i/nb-b/2.d0)*1.d3,(h*j/nh-h/2.d0)*1.d3,0.d0 
      write(12,*) (b*i/nb-b/2.d0)*1.d3,(h*j/nh-h/2.d0)*1.d3,szz(i,j,1) 
      write(13,*) (b*i/nb-b/2.d0)*1.d3,(h*j/nh-h/2.d0)*1.d3,szz(i,j,2) 
      write(14,*) (b*i/nb-b/2.d0)*1.d3,(h*j/nh-h/2.d0)*1.d3,szz(i,j,3) 
   enddo
   write(11,*)
   write(12,*)
   write(13,*)
   write(14,*)
enddo

end
