program plate

! F90/Console application version of plate code
! Limits: 10 material sets, 20 plies

! Equation references are to "Laminated Composite Plates," 
! http://web.mit.edu/course/3/3.11/www/modules/laminates.pdf 

use msimsl	! imsl module

real E_1(10)	! Modulus in fiber direction
real E_2(10)	! Modulus in transverse direction
real G_12(10)	! In-plane shear modulus
real nu_12(10)	! Principal Poisson's ratio
real thk(10)	! Ply thickness

integer k		! Ply index
real theta(20)	! Ply orientation (degrees)	

real eps(6)		! vector of strains and curvatures
real eps0(3)	! Laminate in-plane strains
real kappa(3)	! Laminate curvature
real sig(3)		! Ply stresses in 1-2 frame
real sigbar(3)	! Ply stresses in x-y frame
real z(21)		! Vertical coordinate of bottom of ply

real A(3,3)		! Transformation matrix
real A_inv(3,3)	! Inverse transformation matrix
real C(6,6)		! Laminate compliance matrix
real D_bar(3,3)	! Ply stiffness matrix in x-y frame
real E(6,6)		! Laminate stiffness matrix
real N(6)		! Vector of tractions and moments
real R(3,3)		! Reuter's matrix
real R_inv(3,3)	! Inverse Reuter's matrix
real S(3,3)		! Ply compliance matrix in 1-2 frame
real S_bar(3,3)	! Ply compliance matrix in x-y frame

integer m		! material property set index
integer mat(20)	! ply material index
integer Nplies	! number of plies in layup

data R/1.,3*0.,1.,3*0.,2./
data R_inv/1.,3*0.,1.,3*0.,.5/
data E/36*0./ 

open(unit=10, file='user')

                     
! ----------------------------------------------------------------------					   
! Get material properties

do m = 1,10

  write(10,*) 'Material ', m,': '
  write(10,*) 'Enter modulus in fiber direction (-1 to stop): '
  read(10,*) E_1(m)
    if (E_1(m) .lt. 0.0 ) exit

  write(10,*) 'Enter modulus in transverse direction: '
  read(10,*) E_2(m)

  write(10,*) 'Enter principal Poisson ratio: '
  read(10,*) nu_12(m)

  check = abs((E_1(m)-E_2(m))/E_1(m) )	   ! check for isotropy
  if (check.lt.0.001) then
    G_12(m) = E_1(m)/(2.*(1.+ nu_12(m) ))  ! calculate G_12 if isotropic
  else
    write(10,*) 'Enter shear modulus: '
	read(10,*) G_12(m)  ! otherwise read G_12
  end if				  

  write(10,*) 'Enter ply thickness: '
  read(10,*) thk(m)

end do 

! ----------------------------------------------------------------------					   
! Define layup (ply orientations)

z(1)=0.	  ! vertical coordinate, bottom of first ply 

! loop over plies
do k = 1,20

  write(10,*) 'Define layup sequence, starting at bottom...'
  write(10,*) '   (use negative material set number to stop)'
  write (*,50) k
50  format (/' enter material set number for ply number',i3,': ')
 
  read (10,*) m
    if (m.lt.0) exit
  mat(k) = m

  write(10,*) 'enter ply angle: '
  read(10,*) theta(k)

  z(k+1)=z(k)+thk(m)
  Nplies = k

end do

! Loop over plies again,
! adjust ply coordinates to measure from centerline

z0 = z(Nplies+1)/2.	! midpoint coordinate

do k = 1, Nplies+1
  z(k) = z(k) - z0
end do


! ----------------------------------------------------------------------					   
! Loop over plies, 
! compute contributions to stiffness matrix E(6,6) (Eqn. 25)

do k = 1,Nplies  
  
  z1= z(k+1) - z(k)
  z2= 0.5 * ( z(k+1)**2 - z(k)**2 )
  z3=(1./3.) * ( z(k+1)**3 - z(k)**3)
  call get_D_bar (D_bar)

  do i=1,3
  do j=1,3

    E(i,j)    =  E(i,j) +     D_bar(i,j)*z1	 ! A matrix	(Eqn. 21)
    E(i+3,j)  =  E(i+3,j) +   D_bar(i,j)*z2	 ! B matrix	(Eqn. 22)
    E(i,j+3)  =  E(i,j+3) +   D_bar(i,j)*z2	 ! ", by symmetry
    E(i+3,j+3)=  E(i+3,j+3) + D_bar(i,j)*z3	 ! D matrix	(Eqn. 24)

  end do
  end do

end do 

! ----------------------------------------------------------------------					   
! print stiffness matrix E(6,6)	 

write(10,120)
120 format(/' laminate stiffness matrix:',/)

do i=1,6

  write(10,130) (E(i,j),j=1,6)
130 format (4x,3e12.4,2x,3e12.4)
    if (i.eq.3) write(10,*)

end do

! obtain and print laminate compliance matrix C(6,6)

! call invert (E,6,6,C)	   ! Numerical recipes inversion routine
call LINRG (6,E,6,C,6)	   ! IMSL inversion routine

write(10,310)
310   format(/' laminate compliance matrix:',/)

do i=1,6

  write(10,130) (C(i,j),j=1,6)
  if (i.eq.3) write(10,*)

end do

! ----------------------------------------------------------------------
! obtain traction-moment vector	N(6)

 write(10,*)
 write(10,*) 'input tractions and moments...'
 write(10,*)

 write(10,*) '   Nx: ' ; read(10,*) N(1)
 write(10,*) '   Ny: ' ; read(10,*) N(2)
 write(10,*) '  Nxy: ' ; read(10,*) N(3)
 write(10,*) '   Mx: ' ; read(10,*) N(4)
 write(10,*) '   My: ' ; read(10,*) N(5)
 write(10,*) '  Mxy: ' ; read(10,*) N(6)

! ----------------------------------------------------------------------
! compute and print resulting strains and rotations

eps = matmul (C,N)

write(10,150) (eps(i),i=1,6)
150   format(/' midplane strains:',//3x,'eps-xx =',e12.4,  &
   /3x,'eps-yy =',e12.4,/3x,'eps-xy =',e12.4,			   &
   //' rotations:',//3x,'kappa-xx =',e12.4,				   &
   /3x,'kappa-yy= ',e12.4,/3x,'kappa-xy =',e12.4//)


! ----------------------------------------------------------------------
! compute ply stresses

write(10,160)
160   format (/' stresses:',/2x,'ply',5x,'sigma-1',	 &
        5x,'sigma-2',4x,'sigma-12'/)

do i=1,3  ! get midplane strain and curvature vectors

  eps0(i)  = eps(i)
  kappa(i) = eps(i+3)

end do
  
do k=1,Nplies	! loop over plies (Eqn. 16)

  zctr=(z(k)+z(k+1))/2.		! center coordinate of ply
  call get_D_bar (D_bar)	! this also gets A for current ply

  sigbar = matmul(D_bar,eps0) + zctr * matmul(D_bar,kappa)	! x-y frame
  sig = matmul(A,sigbar)  ! transform to 1-2 frame

  write(10,200) k,sig
  200 format (3x,i2,3e12.4)

end do

write(10,*) 'Enter any number to exit program (screen will vanish): '
read(10,*) k_exit
      
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
!
! internal routines
					   
contains

   subroutine get_A (A)
   ! Transformation matrix for x-y --> 	1-2	(Eqn. 6)
   dimension A(3,3)

   thet = theta(k) * 3.14159/180.

    sc = sin(thet)*cos(thet)
    s2 = (sin(thet))**2
    c2 = (cos(thet))**2

    A(1,1) = c2
    A(2,1) = s2
    A(3,1) = -1.*sc

    A(1,2) = s2
    A(2,2) = c2
    A(3,2) = sc

    A(1,3) = 2.*sc
    A(2,3) = -2.*sc
    A(3,3) = c2 - s2  
	
	end subroutine get_A
	
	! ----------------------------------------------------------------------					   

	subroutine get_A_inv (A_inv)
	! Transformation matrix for 1-2 --> x-y
	dimension A_inv(3,3)

    thet = theta(k) * 3.14159/180.

    sc = sin(thet)*cos(thet)
    s2 = (sin(thet))**2
    c2 = (cos(thet))**2

    A_inv(1,1) = c2
    A_inv(2,1) = s2
    A_inv(3,1) = sc

    A_inv(1,2) = s2
    A_inv(2,2) = c2
    A_inv(3,2) = -1.*sc

    A_inv(1,3) = -2.*sc
    A_inv(2,3) = 2.*sc
    A_inv(3,3) = c2 - s2

	end subroutine get_A_inv

! ----------------------------------------------------------------------					   

  subroutine get_D_bar (D_bar)
  ! Stiffness matrix for current ply in x-y frame
  dimension D_bar(3,3)

  ! First, compute compliance matrix S_bar in x-y frame
  ! [S_bar] = [R][A]-1[R]-1[S][A] (Eqn. 11)
   
  call get_A (A); call get_A_inv (A_inv); call get_S(S)   
  
  S_bar = matmul(R,matmul(A_inv,matmul(R_inv,matmul(S,A))))

  ! invert S_bar to obtain D_bar	

  call LINRG (3,S_bar,3,D_bar,3) 	    ! IMSL inversion routine


  end subroutine get_D_bar

! ----------------------------------------------------------------------

  subroutine get_S (S)
  ! Compliance matrix for current ply in 1-2 frame	(Eqn. 4)
  dimension S(3,3)

	m=mat(k)

    S(1,1) = 1./E_1(m)
    S(2,1) = -nu_12(m) / E_1(m)
    S(3,1) = 0.

    S(1,2) = S(2,1)
    S(2,2) = 1./E_2(m)
    S(3,2) = 0.

    S(1,3) = 0.
    S(2,3) = 0.
    S(3,3) = 1./G_12(m) 
	
   end subroutine get_S

! ----------------------------------------------------------------------					   


end program plate