program BishopHill

!  gfortran -g -O0 -o BishopHill Bishop_Hill-DATE.f90
  !       solve the Bishop Hill model,  adr 28 iv 01
  !  updated for crystal orientation, 21 iv 05
  !  formatted to f90, x 11

  implicit none

  real scalar

  real rtmp

  real stress,strain,dw,dwmax
  integer index,i,j,k,l
  dimension stress(28,6),strain(6)
  real a(3,3),hkl(3),uvw(3),t1(3)
  real eps(3,3),eps_xtal(3,3)
  real strain_xtal(6)
  real rtmp1,rtmp2,tayf
  real evm,enorm
  real sqr3
  real tayf2

  sqr3 = sqrt(3.)
  write(*,*) 'Welcome to the Reid version of the BH algorithm'
  write(*,*) 'Enter the strain as de11,de22,2*de23,2*de31,2*de12'
  read(*,*) strain(1),strain(2),strain(4),strain(5),strain(6)
  strain(3)=0.-strain(1)-strain(2)
  write(*,*) 'Full strain vector: '
  print"(2x,6(2x,f7.3))", (strain(i),i=1,6)

  eps(1,1) = strain(1)
  eps(2,2) = strain(2)
  eps(3,3) = strain(3)
  eps(2,3) = strain(4)/2.
  eps(3,1) = strain(5)/2.
  eps(1,2) = strain(6)/2.
  eps(3,2) = strain(4)/2.
  eps(1,3) = strain(5)/2.
  eps(2,1) = strain(6)/2.
  !	evm = sqrt((2./9.*((strain(1)-strain(2))**2 + 
  !     1     (strain(2)-strain(3))**2 + (strain(3)-strain(1))**2)) +
  !     2     ((eps(1,2)**2 + eps(1,3)**2 + eps(2,3)**2)/3.))
  evm = (2./sqr3)*sqrt(0.5*(eps(1,1)**2 + eps(2,2)**2 + eps(3,3)**2  &
       + eps(2,1)**2 + eps(3,1)**2 + eps(3,2)**2 +  &
       eps(1,2)**2 + eps(1,3)**2 + eps(2,3)**2))
  print"(' eps vM = ',f7.3)",evm

  write(*,*)

  write(*,*) 'Enter (hkl) // sample.3 as 3 numbers'
  read(*,*) hkl(1),hkl(2),hkl(3)
  write(*,"( 'Plane  HKL = ',3(f8.1))") hkl
  call vecnorm(hkl)
10 write(*,*) 'Enter [uvw] // sample.1 as 3 numbers'
  read(*,*) uvw(1),uvw(2),uvw(3)
  write(*,"( 'Plane  UVW = ',3(f8.1))") uvw

  call vecnorm(uvw)
  rtmp=scalar(hkl,uvw)
  if(abs(rtmp).gt.1e-4) then
     write(*,*) 'Sorry, uvw is not perp. to hkl; try again'
     goto 10
  endif
  call vecpro2(hkl,uvw,t1)
  call vecnorm(t1)
  do 20, i=1,3
     a(i,1)=uvw(i)
     a(i,2)=t1(i)
     a(i,3)=hkl(i)
     !	   a(i,3)=t1(i)
     !	   a(i,2)=hkl(i)
20 enddo			!  standard arrangement, not Canova matrix!
  print*,'  Orientation Matrix: '
  do i=1,3
     print"('[ ',3(1x,f7.3),' ]')", (a(i,j),j=1,3)
  enddo

  do i = 1,3
     do j = 1,3
        eps_xtal(i,j) = 0.
        do k =1,3
           do l = 1,3
              eps_xtal(i,j) = eps_xtal(i,j) + a(i,k)*a(j,l)*eps(k,l)
           enddo		!  L
        enddo		!  K
     enddo			!  J
  enddo			!  I

  print*,'  Strain in Xtal coords: '
  do i=1,3
     print"('[ ',3(1x,f7.2),' ]')", (eps_xtal(i,j),j=1,3)
  enddo

  evm = (2./sqr3)*sqrt(0.5*(eps_xtal(1,1)**2 + eps_xtal(2,2)**2  &
       + eps_xtal(3,3)**2 +   &
       eps_xtal(2,1)**2 + eps_xtal(3,1)**2 + eps_xtal(3,2)**2 +  &
       eps_xtal(1,2)**2 + eps_xtal(1,3)**2 + eps_xtal(2,3)**2))
  print"(' eps[xtal] vM = ',f7.3)",evm

  enorm = sqrt((eps_xtal(1,1)**2 + eps_xtal(2,2)**2 + eps_xtal(3,3)**2 +  &
       eps_xtal(2,1)**2 + eps_xtal(3,1)**2 + eps_xtal(3,2)**2 +  &
       eps_xtal(1,2)**2 + eps_xtal(1,3)**2 + eps_xtal(2,3)**2))
  print"(' eps[xtal] vM = ',f7.3)",evm

  strain_xtal(1) = eps_xtal(1,1)
  strain_xtal(2) = eps_xtal(2,2)
  strain_xtal(3) = eps_xtal(3,3)
  strain_xtal(4) = eps_xtal(2,3) + eps_xtal(3,2)
  strain_xtal(5) = eps_xtal(3,1) + eps_xtal(1,3)
  strain_xtal(6) = eps_xtal(1,2) + eps_xtal(2,1)

  print*
  print*,'Strain in crystal axes: '
  print"(2x,6(2x,f7.3))", (strain_xtal(i),i=1,6)

  stress(1,1)=1.
  stress(1,2)=-1.
  stress(1,3)=0.
  stress(1,4)=0.
  stress(1,5)=0.
  stress(1,6)=0.
  !  # 1     
  stress(2,1)=0.
  stress(2,2)=1.
  stress(2,3)=-1.
  stress(2,4)=0.
  stress(2,5)=0.
  stress(2,6)=0.
  !  # 2     
  stress(3,1)=-1.
  stress(3,2)=0.
  stress(3,3)=1.
  stress(3,4)=0.
  stress(3,5)=0.
  stress(3,6)=0.
  !  # 3     
  stress(4,1)=0.
  stress(4,2)=0.
  stress(4,3)=0.
  stress(4,4)=1.
  stress(4,5)=0.
  stress(4,6)=0.
  !  # 4     
  stress(5,1)=0.
  stress(5,2)=0.
  stress(5,3)=0.
  stress(5,4)=0.
  stress(5,5)=1.
  stress(5,6)=0.
  !  # 5     
  stress(6,1)=0.
  stress(6,2)=0.
  stress(6,3)=0.
  stress(6,4)=0.
  stress(6,5)=0.
  stress(6,6)=1.
  !  # 6     
  stress(7,1)=0.5
  stress(7,2)=-1.
  stress(7,3)=0.5
  stress(7,4)=0.
  stress(7,5)=0.5
  stress(7,6)=0.
  !  # 7     
  stress(8,1)=0.5
  stress(8,2)=-1.
  stress(8,3)=0.5
  stress(8,4)=0.
  stress(8,5)=-0.5
  stress(8,6)=0.
  !  # 8     
  stress(9,1)=-1.
  stress(9,2)=0.5
  stress(9,3)=0.5
  stress(9,4)=0.5
  stress(9,5)=0.
  stress(9,6)=0.
  !  # 9     
  stress(10,1)=-1.
  stress(10,2)=0.5
  stress(10,3)=0.5
  stress(10,4)=-0.5
  stress(10,5)=0.
  stress(10,6)=0.
  !  # 10     
  stress(11,1)=0.5
  stress(11,2)=0.5
  stress(11,3)=-1.
  stress(11,4)=0.
  stress(11,5)=0.
  stress(11,6)=0.5
  !  # 11     
  stress(12,1)=0.5
  stress(12,2)=0.5
  stress(12,3)=-1.
  stress(12,4)=0.
  stress(12,5)=0.
  stress(12,6)=-0.5
  !       
  stress(13,1)=0.5
  stress(13,2)=0.
  stress(13,3)=-0.5
  stress(13,4)=0.5
  stress(13,5)=0.
  stress(13,6)=0.5
  !       
  stress(14,1)=0.5
  stress(14,2)=0.
  stress(14,3)=-.5
  stress(14,4)=-0.5
  stress(14,5)=0.
  stress(14,6)=0.5
  !       
  stress(15,1)=0.5
  stress(15,2)=0.
  stress(15,3)=-0.5
  stress(15,4)=0.5
  stress(15,5)=0.
  stress(15,6)=-0.5
  !       
  stress(16,1)=0.5
  stress(16,2)=0.
  stress(16,3)=-0.5
  stress(16,4)=-0.5
  stress(16,5)=0.
  stress(16,6)=-0.5
  !       
  stress(17,1)=0.
  stress(17,2)=-0.5
  stress(17,3)=0.5
  stress(17,4)=0.
  stress(17,5)=0.5
  stress(17,6)=0.5
  !       
  stress(18,1)=0.
  stress(18,2)=-0.5
  stress(18,3)=0.5
  stress(18,4)=0.
  stress(18,5)=-0.5
  stress(18,6)=0.5
  !       
  stress(19,1)=0.
  stress(19,2)=-0.5
  stress(19,3)=0.5
  stress(19,4)=0.
  stress(19,5)=0.5
  stress(19,6)=-0.5
  !       
  stress(20,1)=0.
  stress(20,2)=-0.5
  stress(20,3)=0.5
  stress(20,4)=0.
  stress(20,5)=-0.5
  stress(20,6)=-0.5
  !       
  stress(21,1)=-0.5
  stress(21,2)=0.5
  stress(21,3)=0.
  stress(21,4)=0.5
  stress(21,5)=0.5
  stress(21,6)=0.
  !       
  stress(22,1)=-0.5
  stress(22,2)=0.5
  stress(22,3)=0.
  stress(22,4)=-0.5
  stress(22,5)=0.5
  stress(22,6)=0.
  !       
  stress(23,1)=-0.5
  stress(23,2)=0.5
  stress(23,3)=0.
  stress(23,4)=0.5
  stress(23,5)=-0.5
  stress(23,6)=0.
  !       
  stress(24,1)=-0.5
  stress(24,2)=0.5
  stress(24,3)=0.
  stress(24,4)=-0.5
  stress(24,5)=-0.5
  stress(24,6)=0.
  !       
  stress(25,1)=0.
  stress(25,2)=0.
  stress(25,3)=0.
  stress(25,4)=0.5
  stress(25,5)=0.5
  stress(25,6)=-0.5
  !       
  stress(26,1)=0.
  stress(26,2)=0.
  stress(26,3)=0.
  stress(26,4)=0.5
  stress(26,5)=-0.5
  stress(26,6)=0.5
  !       
  stress(27,1)=0.
  stress(27,2)=0.
  stress(27,3)=0.
  stress(27,4)=-0.5
  stress(27,5)=0.5
  stress(27,6)=0.5
  !       
  stress(28,1)=0.
  stress(28,2)=0.
  stress(28,3)=0.
  stress(28,4)=0.5
  stress(28,5)=0.5
  stress(28,6)=0.5
  !       
  dwmax=0.
  index=0
  do 100, i=1,28
     dw = 0.
     dw = dw - strain_xtal(1)*stress(i,2)
     dw = dw + strain_xtal(2)*stress(i,1)
     dw = dw + strain_xtal(4)*stress(i,4)
     dw = dw + strain_xtal(5)*stress(i,5)
     dw = dw + strain_xtal(6)*stress(i,6)
     !       dW = -Bde11 + Ade22 + 2Fde23 + 2Gde31 +2Hde12
     if(dw.gt.dwmax) then
        dwmax=dw
        index=i
        print*,'dw,i  ',dw,'  ',i
     endif
     if((-1.*dw).gt.dwmax) then
        dwmax=-1.*dw
        index=i+28
        print*,'dw,i  ',dw,'  ',i
     endif
100 end do

  print*
  write(*,*) 'Index of the multiple slip stress state = ',index
  if(index.gt.28) then
     index = index-28
     write(*,*) 'Index of the +/- stress state = ',index
  endif
  print*,'Stress #    A       B       C       F       G       H'
  write(*,120) (stress(index,j),j=1,6)
120 format(6x,6(2x,f6.1))

  rtmp1 = 0.
  do i = 1,6
     rtmp1 = rtmp1 + stress(index,i)**2
  enddo
  !  magnitude of the stress
  rtmp2 = 0.
  do i = 1,6
     rtmp2 = rtmp2 + strain_xtal(i)**2
  enddo
  !  magnitude of the strain in xtal coords.
  !	print*,'dwmax,rtmp1,rtmp2',dwmax,sqrt(rtmp1),sqrt(rtmp2)

  !	tayf = 2.* dwmax / sqrt(rtmp2)
  tayf = dwmax / evm
  !  This is the definition of Taylor factor that allows for multi-axial
  !  strain or stress states, as found in LApp

  print*
  print"('The Taylor factor = ',f7.3,' *sqrt(6)')",tayf
  print"('                  = ',f7.3)",tayf*sqrt(6.)

  !  these next few lines were for curiosity about other normalizations of strain
  !	tayf2 = dwmax / enorm

  !	print*
  !	print"('The Taylor factor = ',f7.3,' *sqrt(6)')",tayf2
  !	print"('                  = ',f7.3)",tayf2*sqrt(6.)

  call exit(0)
end program BishopHill


! ________________________________________

subroutine vecnorm(vec)
 real vec(3),rnorm
 rnorm=0.
 do 10, i=1,3
    rnorm=rnorm+vec(i)**2
10 end do
 if(rnorm.le.0.0) return
 rnorm=sqrt(rnorm)
 do 20, i=1,3
    vec(i)=vec(i)/rnorm
20 end do
 return
end subroutine vecnorm

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

function scalar(a,b)
  !  return SCALAR PRODUCT of A and B
 real scalar,a(3),b(3)
 scalar=0.
 do 100, i=1,3
    scalar=scalar+a(i)*b(i)
100 end do
 return
end function scalar

! ____________________________

subroutine vecpro2(v1,v2,vout)
  !  vector product
 real v1(3),v2(3),vout(3)
 vout(3)=v1(1)*v2(2)-v1(2)*v2(1)
 vout(1)=v1(2)*v2(3)-v1(3)*v2(2)
 vout(2)=v1(3)*v2(1)-v1(1)*v2(3)
 return
end subroutine vecpro2
