c     hotplate - a prompt-driven routine for laminated plate calculations
c                with temperature dependency
 
      real Nx,Ny,Nxy,Mx,My,Mxy
      dimension S(3,3),Sbar(3,3),Dbar(3,3),E(6,7),kwa(6),
     *          T(3,3),Tinv(3,3),R(3,3),Rinv(3,3),Et(6,6),
     *          temp1(3,3),temp2(3,3),temp3(3,3),
     *          eps0(3),xkappa(3),sigbar(3),sig(3),vtemp1(3),
     *          vtemp2(3),E1(30),E2(30),Gnu12(30),G12(30),thk(30),
     *          z(51),mat(50),theta(50),Dsave(3,3,50),Tsave(3,3,50),
     *          cte12(3,30),ctexy(3),ctesave(3,50),
     *          temp4(3),temp5(3),vtemp3(3)
      data R/1.,3*0.,1.,3*0.,2./,Rinv/1.,3*0.,1.,3*0.,.5/,
     *     E/42*0./
 
c----------------------------------------------------------------------
 
      write(6,*)
      write(6,*) 'input temperature difference: '
      read (5,*) deltemp
 
c     input material set selections
 
      i=1
10    write(6,20) i
20    format (' assign properties for lamina type ',i2,'...'/)
 
         write(6,*) 'enter modulus in fiber direction...'
         write(6,*) '   (enter -1 to stop): '
         read (5,*) E1(i)
         if (E1(i) .lt. 0.) go to 30
         write(6,*) 'enter modulus in transverse direction: '
         read (5,*) E2(i)
         write(6,*) 'enter principal Poisson ratio: '
         read (5,*) Gnu12(i)
 
c        check for isotropy
         check=abs((E1(i)-E2(i))/E1(i))
         if (check.lt.0.001) then
            G12(i)=E1(i)/(2.*(1.+Gnu12(i)))
         else
            write(6,*) 'enter shear modulus: '
            read (5,*) G12(i)
         end if
 
         write(6,*) 'enter alpha_1: '
 
         read (5,*) cte12(1,i)
         write(6,*) 'enter alpha_2: '
         read (5,*) cte12(2,i)
         cte12(3,i)=0.
 
         write(6,*) 'enter ply thickness: '
         read (5,*) thk(i)
         i=i+1
 
      go to 10
 
c----------------------------------------------------------------------
c     define layup
 
30    iply=1
      z(1)=0.
 
      write(6,*) 'define layup sequence, starting at bottom...'
      write(6,*) '   (use negative material set number to stop)'
 
40    write (6,50) iply
50    format (/' enter material set number for ply number',i3,': ')
      read (5,*) m
      if (m.lt.0) go to 60
 
         mat(iply)=m
         write(6,*) 'enter ply angle: '
         read (5,*) theta(iply)
         z(iply+1)=z(iply)+thk(m)
         iply=iply+1
 
      go to 40
 
c     compute boundary coordinates (measured from centerline)
 
60    thick=z(iply)
      N = iply-1
      z0 = thick/2.
      np1=N+1
 
      do 70 i=1,np1
         z(i)=z(i)-z0
70    continue
 
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c     loop over plies, form stiffness matrix
 
      do 110 iply=1,N
 
         m=mat(iply)
 
c        form lamina compliance in 1-2 directions (Eqn. 11)
 
 
         S(1,1) = 1./E1(m)
         S(2,1) = -Gnu12(m) / E1(m)
         S(3,1) = 0.
 
         S(1,2) = S(2,1)
         S(2,2) = 1./E2(m)
         S(3,2) = 0.
 
         S(1,3) = 0.
         S(2,3) = 0.
         S(3,3) = 1./G12(m)
 
c----------------------------------------------------------------------
c        transform to x-y axes
c        obtain transformation matrix T (Eqn. 12)
 
 
         thet = theta(iply) * 3.14159/180.
 
         sc = sin(thet)*cos(thet)
         s2 = (sin(thet))**2
         c2 = (cos(thet))**2
 
         T(1,1) = c2
         T(2,1) = s2
         T(3,1) = -1.*sc
 
         T(1,2) = s2
         T(2,2) = c2
         T(3,2) = sc
 
         T(1,3) = 2.*sc
         T(2,3) = -2.*sc
         T(3,3) = c2 - s2
 
c        inverse transformation matrix (Eqn. 13)
 
         Tinv(1,1) = c2
         Tinv(2,1) = s2
         Tinv(3,1) = sc
 
         Tinv(1,2) = s2
         Tinv(2,2) = c2
         Tinv(3,2) = -1.*sc
 
         Tinv(1,3) = -2.*sc
         Tinv(2,3) = 2.*sc
         Tinv(3,3) = c2 - s2
 
c        transformation [Sbar] = [R][T]-1[R]-1[S][T] (Eqn. 16)
 
         call matmul (3,3,3,3,3,3,R,Tinv,temp1)
         call matmul (3,3,3,3,3,3,temp1,Rinv,temp2)
         call matmul (3,3,3,3,3,3,temp2,S,temp3)
         call matmul (3,3,3,3,3,3,temp3,T,Sbar)
 
c        transform thermal expansion vector:
c        alpha_xy =  [R][T]-1[R]-1 alpha_12 = temp2 x cte12
 
         temp4(1)=cte12(1,m)
         temp4(2)=cte12(2,m)
         temp4(3)=cte12(3,m)
         call matmul(3,3,3,3,3,1,temp2,temp4,ctexy)
         ctesave(1,iply)=ctexy(1)
         ctesave(2,iply)=ctexy(2)
         ctesave(3,iply)=ctexy(3)
 
c----------------------------------------------------------------------
c        invert Sbar (transformed compliance matrix) to obtain
c               Dbar (transformed stiffness matrix)
c        start by setting Dbar = Sbar, then call inversion routine
 
         do 80 i=1,3
         do 80 j=1,3
            Dbar(i,j)=Sbar(i,j)
80       continue
 
         call matinv(isol,idsol,3,3,Dbar,3,kwa,det)
 
c        save Dbar and Tinv matrices
 
         do 90 i=1,3
         do 90 j=1,3
            Dsave(i,j,iply)=Dbar(i,j)
            Tsave(i,j,iply)=Tinv(i,j)
90       continue
 
c        add to laminate stiffness matrix
 
         ip1=iply+1
         z1=        (z(ip1)   -z(iply)   )
         z2=    0.5*(z(ip1)**2-z(iply)**2)
         z3=(1./3.)*(z(ip1)**3-z(iply)**3)
         do 100 i=1,3
         do 100 j=1,3
            E(i,j)    =  E(i,j) +     Dbar(i,j)*z1
            xx        =               Dbar(i,j)*z2
            E(i+3,j)  =  E(i+3,j) +   xx
            E(i,j+3)  =  E(i,j+3) +   xx
            E(i+3,j+3)=  E(i+3,j+3) + Dbar(i,j)*z3
100      continue
 
c        compute thermal contribution to load vector
 
         call matmul(3,3,3,3,3,1,Dbar,ctexy,temp5)
         do 101 i=1,3
            E(i  ,7) = E(i  ,7) + temp5(i)*deltemp*z1
            E(i+3,7) = E(i+3,7) + temp5(i)*deltemp*z2
 101     continue
 
c     end loop over plies; stiffness matrix now formed
110   continue
 
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c     output stiffness matrix
 
      write(6,120)
120   format(/' laminate stiffness matrix:',/)
      do 140 i=1,6
         write(6,130) (e(i,j),j=1,6)
130      format (4x,3e12.4,2x,3d12.4)
         if (i.eq.3) write(6,*)
140   continue
 
c----------------------------------------------------------------------
c     obtain and print laminate compliance matrix
 
      do 300 i=1,6
      do 300 j=1,6
         Et(i,j)=E(i,j)
300   continue
 
      call matinv(isol,idsol,6,6,Et,6,kwa,det)
 
      write(6,310)
310   format(/' laminate compliance matrix:',/)
      do 320 i=1,6
        write(6,130) (Et(i,j),j=1,6)
        if (i.eq.3) write(6,*)
320   continue
 
c----------------------------------------------------------------------
c     obtain traction-moment vector
 
      write(6,*)
      write(6,*) 'input tractions and moments...'
      write(6,*)
      write(6,*) '   Nx: '
      read (5,*) Nx
      write(6,*) '   Ny: '
      read (5,*) Ny
      write(6,*) '  Nxy: '
      read (5,*) Nxy
      write(6,*) '   Mx: '
      read (5,*) Mx
      write(6,*) '   My: '
      read (5,*) My
      write(6,*) '  Mxy: '
      read (5,*) Mxy
 
c     add thermal loading
 
      E(1,7) = Nx  + E(1,7)
      E(2,7) = Ny  + E(2,7)
      E(3,7) = Nxy + E(3,7)
      E(4,7) = Mx  + E(4,7)
      E(5,7) = My  + E(5,7)
      E(6,7) = Mxy + E(6,7)
 
 
c----------------------------------------------------------------------
c     solve resulting system; print strains and rotations
 
      call matinv(isol,idsol,6,7,e,6,kwa,det)
      write(6,150) (e(i,7),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//)
 
c----------------------------------------------------------------------
c     compute ply stresses
 
      write(6,160)
160   format (/' stresses:',/2x,'ply',5x,'sigma-1',
     *        5x,'sigma-2',4x,'sigma-12'/)
 
      do 210 iply=1,N
 
         do 180 i=1,3
            eps0(i)  =e(i  ,7)
            xkappa(i)=e(i+3,7)
            ctexy(i) =ctesave(i,iply)
            do 170 j=1,3
               Dbar(i,j)=Dsave(i,j,iply)
               Tinv(i,j)=Tsave(i,j,iply)
170         continue
180      continue
 
         call matmul (3,3,3,3,3,1,Dbar,eps0  ,vtemp1)
         call matmul (3,3,3,3,3,1,Dbar,xkappa,vtemp2)
         call matmul (3,3,3,3,3,1,Dbar,ctexy ,vtemp3)
 
         zctr=(z(iply)+z(iply+1))/2.
         do 190 i=1,3
            sigbar(i) = vtemp1(i) + zctr*vtemp2(i) - deltemp*vtemp3(i)
190      continue
c
         call matmul (3,3,3,3,3,1,T,sigbar,sig)
         write(6,200) iply,sig
200      format (3x,i2,3e12.4)
 
210   continue
 
      stop
      end
 
 
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c --------  library routines for matrix operations -----------
 
      subroutine matmul(lra,lrb,lrc,i,j,k,a,b,c)
c
c     this subroutine performs the multiplication of two
c     two-dimensional matrices (a(i,j)*b(j,k) = c(i,k)).
c
 
c     lra - row dimension of "a" (multiplier) matrix
c     lrb - row dimension of "b" (multiplicand) matrix
c     lrc - row dimension of "c" (product) matrix
c     i   - actual number of rows in "a"
c     j   - actual number of columns in "a"
c     k   - actual number of columns in "b"
c     a   - multiplier
c     b   - multiplicand
c     c   - product
 
      dimension a(1), b(1), c(1)
      do 20 l = 1,i
      nm1 = 0
      lm = l
      do 20 m = 1,k
      c(lm) = 0.0
      nm = nm1
      ln = l
      do 10 n = 1,j
      nm = nm + 1
      c(lm) = c(lm) + a(ln)*b(nm)
 10   ln = ln + lra
      nm1 = nm1 + lrb
 20   lm = lm + lrc
      return
      end
 
 
 
      subroutine matinv(isol,idsol,nr,nc,a,mra,kwa,det)
c
c     this subroutine finds the inverse and/or solves
c     simultaneous equations, or neither, and
c     calculates a determinant of a real matrix.
c
 
c     isol - communications flag (output)
c        1 - inverse found or equations solved
c        2 - unable to solve
c        3 - input error
c     idsol - determinant calculation flag (output)
c        1 - did not overflow
c        2 - overflow
c     nr - number of rows in input matrix "a"
c     nc - number of columns in "a"
c     a  - input matrix, first "nr" columns will be inverted
c             on output, "a" is converted to a-inverse
c     mra - row dimension of "a" matrix
c     kwa - work array
c     det - value of determinant (if idsol = 1)
 
      dimension a(1), kwa(1)
      ir = nr
      isol = 1
      idsol = 1
      if(nr.le.0) go to 330
      if((ir-mra).gt.0) go to 330
      ic = iabs(nc)
      if ((ic - ir).lt.0) ic = ir
      ibmp = 1
      jbmp = mra
      kbmp = jbmp + ibmp
      nes = ir*jbmp
      net = ic*jbmp
      if(nc) 10,330,20
 10   mdiv = jbmp + 1
      iric = ir - ic
      go to 30
 20   mdiv = 1
 30   mad = mdiv
      mser = 1
      kser = ir
      mz = 1
      det = 1.0
 40   piv = 0.
      i = mser
 50   if (( i - kser).gt.0) go to 70
      if((abs(a(i))-piv).le.0.) go to 60
      piv = abs(a(i))
      ip = i
 60   i = i + ibmp
      go to 50
 70   if(piv.eq.0.) go to 340
      if(nc.lt.0) go to 80
      i = ip-((ip - 1)/jbmp)*jbmp
      j = mser - ((mser - 1)/jbmp)*jbmp
      jj = mser/kbmp + 1
      ii = jj + (ip -mser)
      kwa(jj) = ii
      go to 90
 80   i = ip
      j = mser
 90   if (ip - mser) 330,120,100
 100  if ((j - net).gt.0) go to 110
      psto = a(i)
      a(i) = a(j)
      a(j) = psto
      i = i + jbmp
      j = j + jbmp
      go to 100
 110  det = - det
 120  psto = a(mser)
      det = det*psto
 130  if (det.eq.0.) goto 150
 140  psto = 1./psto
      go to 160
 150  idsol = 1
      isol = 2
      return
 160  continue
      a(mser) = 1.0
      i = mdiv
 170  if((i - net).gt.0) go to 180
      a(i) = a(i)*psto
      i = i + jbmp
      go to 170
 180  if((mz - kser).gt.0) go to 210
      if((mz-mser).eq.0) go to 200
      i = mad
      j = mdiv
      psto = a(mz)
      if(psto.eq.0.) go to 200
      a(mz) = 0.
 190  if((j-net).gt.0) go to 200
      a(i) = a(i) - a(j)*psto
      j = j + jbmp
      i = i + jbmp
      go to 190
 200  mad = mad + ibmp
      mz = mz + ibmp
      go to 180
  210 continue
c 210  need a test here.....call overfl(ivf)
c      go to (350,220),ivf
ccccccc  need at test here, anyhow
 220  kser = kser + jbmp
      if ((kser-nes).gt.0) go to 260
      mser = mser + kbmp
      if(nc.lt.0) go to 230
      mdiv = mdiv + ibmp
      mz = ((mser - 1)/jbmp)*jbmp + 1
      mad = 1
      go to 40
 230  mdiv = mdiv + kbmp
      if(iric.ne.0) go to 240
      mz = mser + ibmp
      go to 250
 240  mz = ((mser - 1)/jbmp)*jbmp + 1
 250  mad = mz + jbmp
      go to 40
 260  if(nc.lt.0) return
      jr = ir
 270  if(jr) 330,360,280
 280  if(kwa(jr) - jr) 330,320,290
 290  k = (jr - 1)*jbmp
      j = k + ir
      l = (kwa(jr) - 1)*jbmp + ir
 300  if(j - k) 330,320,310
 310  psto = a(l)
      a(l) = a(j)
      a(j) = psto
      j = j - ibmp
      l = l - ibmp
      go to 300
 320  jr = jr - 1
      go to 270
 330  isol = 3
      return
 340  det = 0.
      isol = 2
      idsol = 1
      return
 350  isol = 2
      idsol = 2
 360  return
      end