inversewavspd.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	inversewavspd (xi, et, c, rho, a)

Function/Subroutine Documentation

inversewavspd()

subroutine inversewavspd	(	real*8, intent(in)	xi,
		real*8, intent(in)	et,
		real*8, dimension(3,3,3,3), intent(in)	c,
		real*8, intent(in)	rho,
		real*8, intent(inout)	a
	)

!
!     Calculates the propagation wave speed in a material, up to its 21 
!     constants. Subroutine for calcmatwavsps.f

!     Based on the procedure in:
!     C. Lane. The Development of a 2D Ultrasonic Array Inspection 
!     for Single Crystal Turbine Blades.
!     Switzerland: Springer International Publishing, 2014.
!
!      CARLO MONJARAZ TEC (CMT)
!
!       INPUT:
!       
!       xi,et: values within a square domain between -1 and 1.
!              correlate in a unique way with 0<=phi<=pi and
!              -pi<=theta<=pi
!
!       elas: c(3,3,3,3) - The elasticity vector
!             
!       rho: double - Density of the material
!        
!        
!       OUTPUT:
!       
!       a: inverse of the wave speed 
!
      implicit none
!     
      integer i,j,k,l,ier,matz,ndim
!
      real*8 c(3,3,3,3),rho,xn(3),cm(3,3,3),a,xi,et,
     &       cmm(3,3),dd,al(3),alz(3,3),fv1(3),fv2(3),
     &       theta,phi,pi,p3(3),v(3),
     &       speed
!
      intent(in) xi,et,c,rho
!
      intent(inout) a
!     
      pi=4.d0*datan(1.d0)
!
      theta=xi*pi
      phi=(et+1.d0)*pi/2.d0
!     
      do l=1,3
         v(l)=0.
         do k=1,3
            cmm(k,l)=0.
            do j=1,3
               cm(j,l,k)=0.
            enddo
         enddo
      enddo            
!     
      xn(1)=dcos(theta)*dsin(phi)
      xn(2)=dsin(theta)*dsin(phi)
      xn(3)=dcos(phi)
!     
!     c ------------ PER EAcH DIREcTION find wave speed-----------------------
!     
      do l=1,3
         do k=1,3
            do i=1,3
               do j=1,3
                  cm(l,k,i)=cm(l,k,i)+c(l,k,j,i)*xn(j)
               enddo
            enddo        
         enddo
      enddo
!     
      do k=1,3
         do i=1,3
            do l=1,3
               cmm(k,i)=cmm(k,i)+cm(l,k,i)*xn(l)
            enddo
         enddo        
      enddo
!     
      ndim=3
      matz=1
      ier=0
!
!     ---------reset vars for EIGvALUES
!
      do j=1,3
         al(j)=0.
         fv1(j)=0.
         fv2(j)=0.
         do i=1,3
            alz(j,i)=0.
         enddo
      enddo
!     
      call rs(ndim,ndim,cmm,al,matz,alz,fv1,fv2,ier)
!     
!           ------normalizing eigenvectors to P vectors----------
!     
      dd=dsqrt(alz(1,3)**2+alz(2,3)**2+alz(3,3)**2)
      p3(1)=alz(1,3)/dd
      p3(2)=alz(2,3)/dd
      p3(3)=alz(3,3)/dd
!     
      do l=1,3
         do k=1,3
            cmm(k,l)=0.
            do j=1,3
               cm(j,l,k)=0.
            enddo
         enddo
      enddo
!     
      do l=1,3
         do j=1,3
            do i=1,3
               do k=1,3
                  cm(l,j,i)=cm(l,j,i)+c(l,k,j,i)*xn(k);
               enddo
            enddo        
         enddo
      enddo
!     
      do  l=1,3
         do  j=1,3        
            do  i=1,3
               cmm(l,j)=cmm(l,j)+cm(l,j,i)*p3(i);        
            enddo
         enddo        
      enddo
!     
      do j=1,3
         do l=1,3
            v(j)=v(j)+cmm(l,j)*p3(l)
         enddo
      enddo
!     
      dd=dsqrt(v(1)**2+v(2)**2+v(3)**2) 
      speed=dd/dsqrt(rho*al(3))
c      speed=dsqrt(dd/rho) 
!     
!     inverse wave speed
!     
      a=1.d0/speed
!     
      return