C squeeze - Compute the forces squeezing on a pair of wires
C           sandwiched between two metal plates.  Each plate
C           is of thickness T and width D.  The two plates
C           are held in contact by screws a distance D apart.
C           Sandwiched between the plates are cylindrical wires
C           of radius R with centers at x1 and x2 relative to
C           the midpoint between the axes of the the two screws.
C           The objective is to compute the forces on the two
C           wires.  As a by-product, the minimum tensions on the
C           two screws required to maintain contact between the
C           plates at +D/2 and -D/2 are also calculated.
C
C author:   Richard Jones, University of Connecticut
C           April 23, 2009
C
C Info:     For information on the properties, see http://www.matweb.com
C           and look at the different aluminum alloys.

      subroutine squeeze(Wcm,Tcm,Dcm,Rcm,x1cm,x2cm)
      real Wcm,Tcm,Dcm,Rcm,x1cm,x2cm   ! all distances in cm
      real W      ! width of the plates (m)
      real T      ! thickness of each plate (m)
      real D      ! distance between the screws holding the plates together (m)
      real R      ! radius of the two wires (m)
      real x1,x2  ! locations of the two wires relative to center of plates (m)
      double precision coef(12)
      common /squeezed/W,T,D,R,x1,x2,coef
      data W,T,D,R,x1,x2/1.0e-2,5.0e-3,1.9e-2,1.0e-5,-2.0e-3,+2.0e-3/
      double precision constrain(12,12)
      integer work(12)
      integer ifail
      integer i,j
      do i=1,12
        do j=1,12
          constrain(i,j) = 0
        enddo
      enddo

      if (Wcm.gt.0) then
        W = Wcm/100
      endif
      if (Tcm.gt.0) then
        T = Tcm/100
      endif
      if (Dcm.gt.0) then
        D = Dcm/100
      endif
      if (Rcm.gt.0) then
        R = Rcm/100
      endif
      if (x1cm.ne.0) then
        x1 = x1cm/100
      endif
      if (x2cm.ne.0) then
        x2 = x2cm/100
      endif
      print *, 'Wire squeeze calculation with wire diameter',
     +         2*R*1e6,' microns'
      print *, 'Plate thickness',T*1e3,' mm, width',W*1e3,' mm,'
      print *, 'and screw spacing',D*1e3,' mm.'

c  1. y=0 at x=D/2
      constrain(1,9) = 1
      constrain(1,10) = D/2
      constrain(1,11) = (D/2)**2
      constrain(1,12) = (D/2)**3
      coef(1) = 0

c  2. y=0 at x=-D/2
      constrain(2,1) = 1
      constrain(2,2) = -D/2
      constrain(2,3) = (D/2)**2
      constrain(2,4) = -(D/2)**3
      coef(2) = 0

c  3. y'=0 at x=D/2
      constrain(3,10) = 1
      constrain(3,11) = 2*(D/2)
      constrain(3,12) = 3*(D/2)**2
      coef(3) = 0

c  4. y'=0 at x=-D/2
      constrain(4,2) = 1
      constrain(4,3) = -2*(D/2)
      constrain(4,4) = 3*(D/2)**2
      coef(4) = 0

c  5. y=R at x=x1
      constrain(5,1) = 1
      constrain(5,2) = x1
      constrain(5,3) = x1**2
      constrain(5,4) = x1**3
      coef(5) = R

c  6. y=R at x=x2
      constrain(6,5) = 1
      constrain(6,6) = x2
      constrain(6,7) = x2**2
      constrain(6,8) = x2**3
      coef(6) = R

c  7. y is continuous at x=x1
      constrain(7,1) = 1
      constrain(7,2) = x1
      constrain(7,3) = x1**2
      constrain(7,4) = x1**3
      constrain(7,5) = -1
      constrain(7,6) = -x1
      constrain(7,7) = -x1**2
      constrain(7,8) = -x1**3
      coef(7) = 0

c  8. y is continuous at x=x2
      constrain(8,5) = 1
      constrain(8,6) = x2
      constrain(8,7) = x2**2
      constrain(8,8) = x2**3
      constrain(8,9) = -1
      constrain(8,10) = -x2
      constrain(8,11) = -x2**2
      constrain(8,12) = -x2**3
      coef(8) = 0

c  9. y' is continuous at x=x1
      constrain(9,2) = 1
      constrain(9,3) = 2*x1
      constrain(9,4) = 3*x1**2
      constrain(9,6) = -1
      constrain(9,7) = -2*x1
      constrain(9,8) = -3*x1**2
      coef(9) = 0

c  10. y' is continuous at x=x2
      constrain(10,6) = 1
      constrain(10,7) = 2*x2
      constrain(10,8) = 3*x2**2
      constrain(10,10) = -1
      constrain(10,11) = -2*x2
      constrain(10,12) = -3*x2**2
      coef(10) = 0

c  11. y'' is continuous at x=x1
      constrain(11,3) = 2
      constrain(11,4) = 6*x1
      constrain(11,7) = -2
      constrain(11,8) = -6*x1
      coef(11) = 0

c  12. y'' is continuous at x=x2
      constrain(12,7) = 2
      constrain(12,8) = 6*x2
      constrain(12,11) = -2
      constrain(12,12) = -6*x2
      coef(12) = 0

c solve for the polynomial coefficients
      call DEQN(12,constrain,12,work,ifail,1,coef)
      if (ifail.ne.0) then
        print *, 'DEQN reports that constrain matrix is singular!'
      endif

c If force on wire 1 is negative, don't hold the plate down!
      if (coef(8).lt.coef(4)) then
        do i=1,12
          do j=1,12
            constrain(i,j) = 0
          enddo
        enddo

c  1a. y=0 at x=D/2
        constrain(1,9) = 1
        constrain(1,10) = D/2
        constrain(1,11) = (D/2)**2
        constrain(1,12) = (D/2)**3
        coef(1) = 0

c  2a. y=0 at x=-D/2
        constrain(2,1) = 1
        constrain(2,2) = -D/2
        constrain(2,3) = (D/2)**2
        constrain(2,4) = -(D/2)**3
        coef(2) = 0

c  3a. y'=0 at x=D/2
        constrain(3,10) = 1
        constrain(3,11) = 2*(D/2)
        constrain(3,12) = 3*(D/2)**2
        coef(3) = 0

c  4a. y'=0 at x=-D/2
        constrain(4,2) = 1
        constrain(4,3) = -2*(D/2)
        constrain(4,4) = 3*(D/2)**2
        coef(4) = 0

c  5a. y''' is continuous at x=x1
        constrain(5,4) = 1
        constrain(5,8) = -1
        coef(5) = 0

c  6a. y=R at x=x2
        constrain(6,5) = 1
        constrain(6,6) = x2
        constrain(6,7) = x2**2
        constrain(6,8) = x2**3
        coef(6) = R

c  7a. y is continuous at x=x1
        constrain(7,1) = 1
        constrain(7,2) = x1
        constrain(7,3) = x1**2
        constrain(7,4) = x1**3
        constrain(7,5) = -1
        constrain(7,6) = -x1
        constrain(7,7) = -x1**2
        constrain(7,8) = -x1**3
        coef(7) = 0

c  8a. y is continuous at x=x2
        constrain(8,5) = 1
        constrain(8,6) = x2
        constrain(8,7) = x2**2
        constrain(8,8) = x2**3
        constrain(8,9) = -1
        constrain(8,10) = -x2
        constrain(8,11) = -x2**2
        constrain(8,12) = -x2**3
        coef(8) = 0

c  9a. y' is continuous at x=x1
        constrain(9,2) = 1
        constrain(9,3) = 2*x1
        constrain(9,4) = 3*x1**2
        constrain(9,6) = -1
        constrain(9,7) = -2*x1
        constrain(9,8) = -3*x1**2
        coef(9) = 0

c  10a. y' is continuous at x=x2
        constrain(10,6) = 1
        constrain(10,7) = 2*x2
        constrain(10,8) = 3*x2**2
        constrain(10,10) = -1
        constrain(10,11) = -2*x2
        constrain(10,12) = -3*x2**2
        coef(10) = 0

c  11a. y'' is continuous at x=x1
        constrain(11,3) = 2
        constrain(11,4) = 6*x1
        constrain(11,7) = -2
        constrain(11,8) = -6*x1
        coef(11) = 0

c  12a. y'' is continuous at x=x2
        constrain(12,7) = 2
        constrain(12,8) = 6*x2
        constrain(12,11) = -2
        constrain(12,12) = -6*x2
        coef(12) = 0

c solve for the polynomial coefficients
        call DEQN(12,constrain,12,work,ifail,1,coef)
        if (ifail.ne.0) then
          print *, 'DEQN reports that constrain matrix is singular!'
        endif

      elseif (coef(12).lt.coef(8)) then
        do i=1,12
          do j=1,12
            constrain(i,j) = 0
          enddo
        enddo

c  1b. y=0 at x=D/2
        constrain(1,9) = 1
        constrain(1,10) = D/2
        constrain(1,11) = (D/2)**2
        constrain(1,12) = (D/2)**3
        coef(1) = 0

c  2b. y=0 at x=-D/2
        constrain(2,1) = 1
        constrain(2,2) = -D/2
        constrain(2,3) = (D/2)**2
        constrain(2,4) = -(D/2)**3
        coef(2) = 0

c  3b. y'=0 at x=D/2
        constrain(3,10) = 1
        constrain(3,11) = 2*(D/2)
        constrain(3,12) = 3*(D/2)**2
        coef(3) = 0

c  4b. y'=0 at x=-D/2
        constrain(4,2) = 1
        constrain(4,3) = -2*(D/2)
        constrain(4,4) = 3*(D/2)**2
        coef(4) = 0

c  5b. y''' is continuous at x=x2
        constrain(5,8) = 1
        constrain(5,12) = -1
        coef(5) = 0

c  6b. y=R at x=x1
        constrain(6,1) = 1
        constrain(6,2) = x1
        constrain(6,3) = x1**2
        constrain(6,4) = x1**3
        coef(6) = R

c  7b. y is continuous at x=x1
        constrain(7,1) = 1
        constrain(7,2) = x1
        constrain(7,3) = x1**2
        constrain(7,4) = x1**3
        constrain(7,5) = -1
        constrain(7,6) = -x1
        constrain(7,7) = -x1**2
        constrain(7,8) = -x1**3
        coef(7) = 0

c  8b. y is continuous at x=x2
        constrain(8,5) = 1
        constrain(8,6) = x2
        constrain(8,7) = x2**2
        constrain(8,8) = x2**3
        constrain(8,9) = -1
        constrain(8,10) = -x2
        constrain(8,11) = -x2**2
        constrain(8,12) = -x2**3
        coef(8) = 0

c  9b. y' is continuous at x=x1
        constrain(9,2) = 1
        constrain(9,3) = 2*x1
        constrain(9,4) = 3*x1**2
        constrain(9,6) = -1
        constrain(9,7) = -2*x1
        constrain(9,8) = -3*x1**2
        coef(9) = 0

c  10b. y' is continuous at x=x2
        constrain(10,6) = 1
        constrain(10,7) = 2*x2
        constrain(10,8) = 3*x2**2
        constrain(10,10) = -1
        constrain(10,11) = -2*x2
        constrain(10,12) = -3*x2**2
        coef(10) = 0

c  11b. y'' is continuous at x=x1
        constrain(11,3) = 2
        constrain(11,4) = 6*x1
        constrain(11,7) = -2
        constrain(11,8) = -6*x1
        coef(11) = 0

c  12b. y'' is continuous at x=x2
        constrain(12,7) = 2
        constrain(12,8) = 6*x2
        constrain(12,11) = -2
        constrain(12,12) = -6*x2
        coef(12) = 0

c solve for the polynomial coefficients
        call DEQN(12,constrain,12,work,ifail,1,coef)
        if (ifail.ne.0) then
          print *, 'DEQN reports that constrain matrix is singular!'
        endif
      endif
      end

      subroutine report
      real W,T,D,R,x1,x2
      double precision coef(12)
      common /squeezed/W,T,D,R,x1,x2,coef
      real B               ! Young's modulus of the plate materials (Pascals)
      parameter (B=69e9)   ! for aluminum
      real A,S1,S2,S3,S4
      A = 2*B*W*T**3
      S1 = A*coef(4)
      S2 = A*(coef(8)-coef(4))
      S3 = A*(coef(12)-coef(8))
      S4 = -A*coef(12)
      print *, 'Forces on wire 1,2 are ',S2,',',S3,' in Newtons'
      print *, 'Forces on screws 1,2 are ',S1,',',S4,' in Newtons'
      end

      real function ycm(xcm) ! (in cm)
      real xcm
      real W,T,D,R,x1,x2
      double precision coef(12)
      common /squeezed/W,T,D,R,x1,x2,coef
      real x,y
      x = xcm/100
      if (abs(x).gt.D/2) then
        y = 0
      elseif (x.lt.x1) then
        y = coef(1)+coef(2)*x+coef(3)*x**2+coef(4)*x**3 
      elseif (x.lt.x2) then
        y = coef(5)+coef(6)*x+coef(7)*x**2+coef(8)*x**3 
      else
        y = coef(9)+coef(10)*x+coef(11)*x**2+coef(12)*x**3 
      endif
      ycm = y*100
      end

#if PROGRAM
      program main
      call squeeze(0,0,0,0,0,0)
      end
#endif