c       case2_step.for

c       calculates stress and deflection in a stepped endplate 

c       Roark & Young, table 24, case 2, with a step

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      dimension amat(4,4), cvec(4), vbeg(4), oswire(3,99), area(99),
     .          twire(3), Ewire(3), dwire(3), EAwire(3), fwire(3),
     .          loadfac(3,99), ncelltot(99), nwire(3), dhole(3),
     .          dzf(2,99)

c          bi = inner radius of region i
      data bi / 0.235, ! inner radius of plate
     .          0.2530, 0.2662, 0.2781, 0.2899, 0.3016, ! superlayer 1
     .          0.3124, 0.3243, 0.3361, 0.3479, 0.3596, ! superlayer 2
     .          0.3646, 0.3764, 0.3882, 0.4000, 0.4118, ! superlayer 3
     .          0.4168, 0.4286, 0.4404, 0.4522, 0.4639, ! superlayer 4
     .          0.4747, 0.4867, 0.4984, 0.5101, 0.5219, ! superlayer 5
     .          0.5269, 0.5387, 0.5505, 0.5623, 0.5741, ! superlayer 6
     .          0.5791, 0.5909, 0.6027, 0.6145, 0.6262, ! superlayer 7
     .          0.6370, 0.6489, 0.6606, 0.6724, 0.6842, ! superlayer 8
     .          0.6892, 0.7010, 0.7128, 0.7246, 0.7364, ! superlayer 9
     .          0.7414, 0.7532, 0.7650, 0.7768, 0.7900, ! superlayer 10
     .          0.8090, 47*0.0 / ! outer radius of plate

c           ihole = 1 -- region has holes
      data ihole / 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1,
     .             0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 49*0 /

c         layer = drift chamber layer number
      data layer / 0,  1,  2,  3,  4, 0,  5,  6,  7,  8, 
     .             0,  9, 10, 11, 12, 0, 13, 14, 15, 16, 
     .             0, 17, 18, 19, 20, 0, 21, 22, 23, 24, 
     .             0, 25, 26, 27, 28, 0, 29, 30, 31, 32, 
     .             0, 33, 34, 35, 36, 0, 37, 38, 39, 40,
     .             49*0 /

c         desired wire tensions in grams (sense, clear, field)
      data twire / 34., 86., 182. /  ! nominal
c      data twire / 34., 86., 86. /   ! if use 80 micron wires for field wires
c       data twire / 3*0.0 /  ! no load case

c         modulus of wire (Pa)
      data Ewire / 3.55e11, 2*7.0e10 /

c         diameter of wires (m)
      data dwire / 0.00002, 0.00008, 0.00012 /

c         diameter of holes (m)
      data dhole / 0.0045, 2*0.0025 /

c         cells per layer
      data ncelltot / 0, 4*96, 0, 4*112, 0, 4*128, 0, 4*144, 0, 4*176,
     .             0, 4*192, 0, 4*208, 0, 4*224, 0, 4*240, 0, 4*256,
     .             49*0 /

c         number of wires per cell (sense, clear, field)
      data loadfac /
     .              3*0, 1,2,3, 1,0,2, 1,0,2, 1,2,2, ! 1 extra field wire 
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, ! in first and last layer
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2,
     .              3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,3,
     .              147*0/

      pi = 4.*atan(1.)

c       set up wire parameters
      do i = 1,3
        fwire(i) = twire(i)*9.8/1000. ! tension in N
        EAwire(i) = Ewire(i)*pi*dwire(i)**2/4.
        nwire(i) = 0
      enddo ! i

      zlen = 2.764 ! wire length (m)

c      setup for plate divided into 51 regions

      n = 51

      rstep = 0.47 ! radius of step in front endplate
      
      t1r = 0.024 ! thickness of rear endplate
      t1f = 0.024 ! thickness for r < rstep
c      t1f = 0.012
      t2f = 0.012 ! thickness for r > rstep
c      t2f = 0.024

c      delta = z offset between inner and outer radius, + means outwards
      delta = 0.0031
c      delta = 0.

      Eal = 7.0e10 ! Young's modulus for 6061aluminum
c      Eal = 7.7e10 ! Young's modulus for 8090 aluminum
      anu = 0.33 ! Poissons' ratio
c        values of holefac and stresscon chosen to match Hodgson's results 
      holefac = 1./1.5 ! modulus reduction factor for holes
c      holefac = 1.
      stressfac = 4.1 ! stress concentration factor around holes
c      stressfac = 1.

      aload = 0.
      do i = 1,51
        rmid(i) = 0.5*(bi(i+1) + bi(i))
        area(i) = pi*(bi(i+1)**2 - bi(i)**2)
        E(i) = Eal
        anui(i) = anu
        t(1,i) = t1f ! front endplate
        if (bi(i) .gt. rstep) t(1,i) = t2f
        t(2,i) = t1r ! rear endplate
        qi(i) = 0.
        if (ihole(i). eq. 1) then
c          holes(i) = holefac
          stresscon(i) = stressfac
          awire = 0.
          do j = 1,3
            loadfac(j,i) = loadfac(j,i)*ncelltot(i) ! convert to wires per layer
            ! calculate area of wires in this layer
            awire = awire + loadfac(j,i)*pi*dhole(j)**2/4.
            nwire(j) = nwire(j) + loadfac(j,i)
            oswire(j,i) = twire(j) ! oswire stores wire tension in grams
            temp = oswire(j,i)*loadfac(j,i)/1000. ! wire tension in N
            aload = aload + temp ! load in kg
            qi(i) = qi(i) + temp*9.8/area(i) ! pressure in Pa


c  turn on the next line for no load conditions
c            qi(i) = 0.


c            write(1,9898) j, i, loadfac(j,i), oswire(j,i), temp
9898         format (' j, i',2i3,' loadfac, oswire, temp ',i5,2e11.3)
          enddo ! j
          ! set hole factor to fraction of area not in holes
c          holes(i) = 1. - awire/area(i)      
          ! set hole factor to holefac
          holes(i) = holefac
        else
            holes(i) = 1.
            stresscon(i) = 1.
        endif

      enddo ! i

      write (1,1) zlen, n, delta
    1 format (' Stress analysis of a chamber ',f5.3,' m long, with ',
     . i3,' regions'/' z offset between inner and outer radii = ',
     . f7.4/)
     
      write (1,2) aload, nwire
    2 format (' wire load = ',f5.0,' kg for ', i5,' sense, ',i5,
     . ' clearing, and ',i5,' field wires'//
     . ' Parameters for the regions:'//
     . ' region, b(cm), P(kPa), E(Gpa), nu, holefac, stress tf(mm),'
     . '  Df,  tr,  Dr')

      do i = 1,n

        do j = 1,2
        Di(j,i) = E(i)*holes(i)*t(j,i)**3/12./(1. - anui(i)**2) ! plate constant
        enddo ! j

        write (1,11) i, bi(i)*100., qi(i)/1000., 
     .  E(i)*1.e-9, anui(i), holes(i), stresscon(i), 
     .  t(1,i)*1000., Di(1,i), t(2,i)*1000., Di(2,i)
   11   format (i4,f10.1,f8.2,f8.2,f6.2,f7.3,f8.2,f7.1,f7.0,f7.1,f7.0)

c         write on unit 2 for TeX table
        write (2,112)  i, bi(i)*100., (loadfac(j,i),j=1,3),
     .    qi(i)/1000., holes(i), stresscon(i), t(1,i)*1000., Di(1,i)
  112   format (i3,' & ',f4.1,' & ',i3,' & ',i3,' & ',i3,' & ',
     .     f4.1,' & ',f4.2,' & ',f3.1,' & ',f3.0,' & ',f7.0,' \t')

   10   format (
     .  ' region',t30,i3/
     .  ' inner radius',t30,1pg11.3,' m'/
     .  ' thickness, b < r < c',t30,g11.3,' m'/
     .  ' pressure',t30,g11.3,' Pascals'/
     .  ' Youngs modulus',t30,1pg11.3,' Pascals'/
     .  ' Poissons ratio',t30,g11.3/
     .  ' Modulus loss to holes',t30,g11.3/
     .  ' Stress conc. at holes',t30,g11.3/
     .  ' D = plate constant',t30,g11.3/)

      enddo ! i
      
      write (1,12) bi(n+1)
   12 format (
     . ' outer radius',t30,g11.3,' m'/)

c      front endplate

c        build the transformation from inner to outer regions
      call build(amat,cvec,1)

c        solve for values at inner radius, case 2c
      call case2c(vbeg,amat,cvec)

c        print final values, and store final displacements in dzf
      call output(1,vbeg,0,1,dzf)

c         case2d, clamped at inner radius, simply supported at outer
c      call case2d(vbeg,amat,cvec)
c      call output(2,vbeg,0,1,dzf)

c      stop
9999  continue

c      rear endplate
      call build(amat,cvec,2)
      call case2c(vbeg,amat,cvec)
      call output(1,vbeg,0,2,dzf)

c      call case2d(vbeg,amat,cvec)
c      call output(2,vbeg,0,2,dzf)


      stop
      end

      subroutine build(amat,cvec,ifr)

c        build transformation between inner and outer radii

c      vec_j(1) -- y at inner edge of region j
c      vec_j(2) -- theta
c      vec_j(3) -- radial bending moment, M_r
c      vec_j(4) -- shear force density, Q

c         form of transformation of vec_j from the inner edge of region j to 
c         vec_i at the outer edge of region j:
c      vec_i(k) = sum_l amat_j(k,l)*vec_j(l) + cvec_j(k),
c          where amat_j and cvec_j are obtained from routine smatrix

c         we also want to transform from the inner edge of region 1 to the
c         outer edge of region j:
c      vec_i(k) = sum_l amat_i(k,l)*vec_j(l) + cvec_i(k)

c         therefore
c      cvec_i(k) = sum_l amat_j(k,l)*cvec_i-1(l) + cvec_j(k),
c      amat_i(k,l) = sum_m amat_j(k,m)*amat_i-1(m,l).

c         I store cvec_i-1 in cvec, and amat_i-1 in amat

c       ifr = 1 -- front, = 2 -- rear

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      dimension amat(4,4), cvec(4), amatj(4,4), cvecj(4), amati(4,4),
     .          cveci(4)


      do k = 1,4
        cvec(k) = 0.
        do l = 1,4
          amat(k,l) = 0.
          if (k .eq. l) amat(k,l) = 1.
        enddo ! l
      enddo ! k

      do j = 1,n

        ! get coeffs at outer edge of region
        call smatrix(j,bi(j+1),amatj,cvecj,ifr) 

c        write (1,93) j, bi(j+1), cvecj, ((amatj(k,l),l=1,4),k=1,4)
   93   format (' region ',i3,' outer radius ',f8.4/
     .   ' cvec',t10,4e11.3/' amat ',(t10,4e11.3))
 
        do k = 1,4
          cveci(k) = cvecj(k) 
          do l = 1,4
            cveci(k) = cveci(k) + amatj(k,l)*cvec(l)
            amati(k,l) = 0.
            do m = 1,4
              amati(k,l) = amati(k,l) + amatj(k,m)*amat(m,l)
            enddo ! m
          enddo ! l
        enddo ! k

        do k = 1,4
          cvec(k) = cveci(k)
          do l = 1,4
            amat(k,l) = amati(k,l)
          enddo ! l
        enddo ! k
      enddo ! j

      return
      end

      subroutine case2c(vbeg,amat,cvec)

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      dimension amat(4,4), cvec(4), vbeg(4), vfin(4)

c       apply B.C.'s at inner and outer radii -- depending on case

c        case 2c, simply supported at both a and b

c        height = 0 at inner and delta atouter radii
      vbeg(1) = 0.
      vfin(1) = delta
c        bend moment = 0 at inner and outer radii
      vbeg(3) = 0.
      vfin(3) = 0.

c        solve for vbeg(2) and veg(4) using
c     vfin(1) = delta = amat(1,2)*vbeg(2) + amat(1,4)*vbeg(4) + cvec(1)
c     vfin(3) = 0 = amat(3,2)*vbeg(2) + amat(3,4)*vbeg(4) + cvec(3)

      denom = amat(1,2)*amat(3,4) - amat(3,2)*amat(1,4)
      vbeg(2) = -((cvec(1)-delta)*amat(3,4) - cvec(3)*amat(1,4))/denom
      vbeg(4) = -(amat(1,2)*cvec(3) - amat(3,2)*(cvec(1)-delta))/denom

      write (1,92) amat(1,2), amat(1,4), cvec(1),
     .             amat(3,2), amat(3,4), cvec(3)
   92 format (' M12 ',e11.3,' M14 ',e11.3,' C1 ',e11.3/
     .        ' M32 ',e11.3,' M34 ',e11.3,' C3 ',e11.3/)

      return
      end

      subroutine case2d(vbeg,amat,cvec)

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      dimension amat(4,4), cvec(4), vbeg(4), vfin(4)

c       apply B.C.'s at inner and outer radii -- depending on case

c        case 2d, fixed at inner radius, simply supported at outer radius

c        height = 0 at inner and outer radii
      vbeg(1) = 0.
      vfin(1) = 0.
c        slope = 0 at inner radius
      vbeg(2) = 0.
c        bend moment = 0 at outer radius
      vfin(3) = 0.

c        solve for vbeg(3) and veg(4) using
c     vfin(1) = 0 = amat(1,3)*vbeg(3) + amat(1,4)*vbeg(4) + cvec(1)
c     vfin(3) = 0 = amat(3,3)*vbeg(3) + amat(3,4)*vbeg(4) + cvec(3)

      denom = amat(1,3)*amat(3,4) - amat(3,3)*amat(1,4)
      vbeg(3) = -(cvec(1)*amat(3,4) - cvec(3)*amat(1,4))/denom
      vbeg(4) = -(amat(1,3)*cvec(3) - amat(3,3)*cvec(1))/denom

      return
      end

      subroutine output(icase,vbeg,iout,ifr,dz)

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      dimension vbeg(4), vbegi(4,99), vfin(4), dz(2,99)
      character cases(2) / 'c', 'd' /
    
      pi = 4.*atan(1.)

      if (iout .eq. 0 .and. ifr. eq. 1) write(1,2)
    2 format (' Results for Front Endplate')
      if (iout .eq. 0 .and. ifr. eq. 2) write(1,3)
    3 format (' Results for Rear Endplate')

      if (iout .eq. 0) write (1,1) cases(icase)
    1 format (/' case 2',a1//
     . ' Values at inner edges of regions' //
     . ' region, y(mm), theta(rad), M_r, Load(kg)')


c        find values at beginning of each region
      do j = 1,4
        vbegi(j,1) = vbeg(j)
      enddo ! j

      do i = 1,n+1
        if (i. gt. 1) 
     .     call values(vbegi(1,i),vbegi(1,i-1),bi(i),i-1,ifr)

c      write (1,91) (vbegi(j,i-1),j=1,4), (vbegi(j,i),j=1,4)
   91 format (' vbeg ',4e11.3/' vfin ',4e11.3)

        alb = 2.*pi*bi(i)*vbegi(4,i)/9.8
        if (iout .eq. 0) write (1,11) i, vbegi(1,i)*1000., vbegi(2,i), 
     .                                vbegi(3,i), alb
   11   format (i4,f10.3,f9.4,f9.1,f8.0)

        enddo ! i

c        print values at middles of regions that have wires

      if (iout .eq. 0) write (1,29)
   29 format (/' Values at middle of regions that have wires'/)

      if (iout. eq. 0) write (1,21)
   21 format(' layer, r(cm), y(mm), theta(rad), M_r(N-m/m), M_t,'
     . ' Q(N), sigma_r(MPa), sigma_t')

      do i = 1,n
          r = rmid(i)
          call values(vfin,vbegi(1,i),r,i,ifr)
c           store displacement; change sign!!!!
          dz(ifr,i) = -vfin(1) 

          amt = vfin(2)*Di(ifr,i)*(1. - anui(i)**2)/r + 
     .          anui(i)*vfin(3)
          sigmar = 6.*vfin(3)/t(ifr,i)**2*1.e-6*stresscon(i) ! include stress 
          sigmat = 6.*amt/t(ifr,i)**2*1.e-6*stresscon(i) ! concentration factor

          if (ihole(i). eq. 1 .and. iout .eq. 0) 
     .       write (1,30) layer(i), r*100., 
     .                 vfin(1)*1000., vfin(2), 
     .                 vfin(3), amt, vfin(4), sigmar, sigmat

c        write on unit 2 for TeX table
        write (2,31) i, r*100., 
     .                 vfin(1)*1000., vfin(2)*1000., 
     .       vfin(3), amt, vfin(4)*2.*pi*bi(i)/9.8, sigmar, sigmat
   31   format (i3,' & ',f5.1,' & ',f6.3,' & ',f6.1,' & ',f7.1,' & ',
     .  f7.1,' & ',f7.0,' & ',f6.1,' & ',f6.1,' \t')

      enddo ! i

      return

9999  continue

c        print values for 100 radii
      do i = 1,101
        r = bi(1) + 0.01*(i- 1)*(bi(n+1)  - bi(1))

        ! find which region we are in
        j = 0
  100   j = j + 1
        if (r .gt. bi(j+1)) go to 100
        if (j .gt. n) j = n
        call values(vfin,vbegi(1,j),r,j,ifr)
        amt = vfin(2)*Di(ifr,j)*(1. - anui(j)**2)/r + anui(j)*vfin(3)
        sigmar = 6.*vfin(3)/t(ifr,j)**2*1.e-6*stresscon(j) ! include stress 
        sigmat = 6.*amt/t(ifr,j)**2*1.e-6*stresscon(j) ! concentration factor
  
        write (1,30) j, r*100., vfin(1)*1000., vfin(2), vfin(3), 
     .   amt, vfin(4), sigmar, sigmat
   30   format (1x,i3, f7.2,',',f8.3,',',f9.4,',',f9.1,',',
     .  f9.1,',',f9.0,',',f9.2,',',f9.2)

      enddo ! i

      return
      end

      subroutine values(vfin,vbeg,r,i,ifr)

      implicit real*8 (a-h), (o-z)

      dimension vfin(4), vbeg(4), amat(4,4), cvec(4)

      call smatrix(i,r,amat,cvec,ifr)

      do j = 1,4

        vfin(j) = cvec(j)

        do k = 1,4
          vfin(j) = vfin(j) + amat(j,k)*vbeg(k)
        enddo !k
      enddo ! j

      return
      end

      subroutine smatrix(i,r,amat,cvec,ifr)

      implicit real*8 (a-h), (o-z)

      common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99),
     .             holes(99), stresscon(99), n, rmid(99), ihole(99), 
     .             layer(99), delta

      common/nu/anu
      dimension amat(4,4), cvec(4)

      anu = anui(i)
      b = bi(i)
      D = Di(ifr,i)
      q = qi(i)

      amat(1,1) = 1.
      amat(1,2) = r*f1(r,b)
      amat(1,3) = r**2/D*f2(r,b)
      amat(1,4) = r**3/D*f3(r,b)
      amat(2,1) = 0.
      amat(2,2) = f4(r,b)
      amat(2,3) = r/D*f5(r,b)
      amat(2,4) = r**2/D*f6(r,b)
      amat(3,1) = 0.
      amat(3,2) = D/r*f7(r,b)
      amat(3,3) = f8(r,b)
      amat(3,4) = r*f9(r,b)
      amat(4,1) = 0.
      amat(4,2) = 0.
      amat(4,3) = 0.
      amat(4,4) = b/r

      cvec(1) = -q*r**4/D*f11(r,b)
      cvec(2) = -q*r**3/D*f14(r,b)
      cvec(3) = -q*r**2*f17(r,b)
      cvec(4) = -q*(r**2 - b**2)/2./r

      return
      end

      real*8 function f1(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f1 = 0.5*(1. + anu)*b/r*log(r/b) + 
     .     0.25*(1. - anu)*(r/b - b/r)
      return
      end

      real*8 function f2(r,b)
      implicit real*8 (a-h), (o-z)
      f2 = 0.25*(1. - (b/r)**2*(1. + 2.*log(r/b)))
      return
      end

      real*8 function f3(r,b)
      implicit real*8 (a-h), (o-z)
      f3 = 0.25*b/r*((1. + (b/r)**2)*log(r/b) + (b/r)**2 - 1.)
      return
      end

      real*8 function f4(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f4 = 0.5*((1. + anu)*b/r + (1. - anu)*r/b)
      return
      end

      real*8 function f5(r,b)
      implicit real*8 (a-h), (o-z)
      f5 = 0.5*(1. - (b/r)**2)
      return
      end

      real*8 function f6(r,b)
      implicit real*8 (a-h), (o-z)
      f6 = 0.25*b/r*((b/r)**2 - 1. + 2.*log(r/b))
      return
      end

      real*8 function f7(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f7 = 0.5*(1. - anu**2)*(r/b - b/r)
      return
      end

      real*8 function f8(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f8 = 0.5*(1. + anu + (1. - anu)*(b/r)**2)
      return
      end

      real*8 function f9(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f9 = b/r*(0.5*(1. + anu)*log(r/b) + 0.25*(1. - anu)*
     .      (1. - (b/r)**2))
      return
      end

      real*8 function f11(r,b)
      implicit real*8 (a-h), (o-z)
      f11 = (1. + 4.*(b/r)**2 - 5.*(b/r)**4 - 4.*(b/r)**2*
     .       (2. + (b/r)**2)*log(r/b))/64.
      return
      end

      real*8 function f14(r,b)
      implicit real*8 (a-h), (o-z)
      f14 = (1. - (b/r)**4 - 4.*(b/r)**2*log(r/b))/16. 
      return
      end

      real*8 function f17(r,b)
      implicit real*8 (a-h), (o-z)
      common/nu/anu
      f17 = 0.25*(1. - 0.25*(1. - anu)*(1. - (b/r)**4) -
     .     (b/r)**2*(1. + (1. + anu)*log(r/b)))
      return
      end