1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine inversewavspd(xi,et,c,rho,a)
20!
21!     Calculates the propagation wave speed in a material, up to its 21
22!     constants. Subroutine for calcmatwavsps.f
23
24!     Based on the procedure in:
25!     C. Lane. The Development of a 2D Ultrasonic Array Inspection
26!     for Single Crystal Turbine Blades.
27!     Switzerland: Springer International Publishing, 2014.
28!
29!      CARLO MONJARAZ TEC (CMT)
30!
31!       INPUT:
32!
33!       xi,et: values within a square domain between -1 and 1.
34!              correlate in a unique way with 0<=phi<=pi and
35!              -pi<=theta<=pi
36!
37!       elas: c(3,3,3,3) - The elasticity vector
38!
39!       rho: double - Density of the material
40!
41!
42!       OUTPUT:
43!
44!       a: inverse of the wave speed
45!
46      implicit none
47!
48      integer i,j,k,l,ier,matz,ndim
49!
50      real*8 c(3,3,3,3),rho,xn(3),cm(3,3,3),a,xi,et,
51     &       cmm(3,3),dd,al(3),alz(3,3),fv1(3),fv2(3),
52     &       theta,phi,pi,p3(3),v(3),
53     &       speed
54!
55!
56!
57      pi=4.d0*datan(1.d0)
58!
59      theta=xi*pi
60      phi=(et+1.d0)*pi/2.d0
61!
62      do l=1,3
63         v(l)=0.d0
64         do k=1,3
65            cmm(k,l)=0.d0
66            do j=1,3
67               cm(j,l,k)=0.d0
68            enddo
69         enddo
70      enddo
71!
72      xn(1)=dcos(theta)*dsin(phi)
73      xn(2)=dsin(theta)*dsin(phi)
74      xn(3)=dcos(phi)
75!
76!     c ------------ PER EAcH DIREcTION find wave speed-----------------------
77!
78      do l=1,3
79         do k=1,3
80            do i=1,3
81               do j=1,3
82                  cm(l,k,i)=cm(l,k,i)+c(l,k,j,i)*xn(j)
83               enddo
84            enddo
85         enddo
86      enddo
87!
88      do k=1,3
89         do i=1,3
90            do l=1,3
91               cmm(k,i)=cmm(k,i)+cm(l,k,i)*xn(l)
92            enddo
93         enddo
94      enddo
95!
96      ndim=3
97      matz=1
98      ier=0
99!
100!     ---------reset vars for EIGvALUES
101!
102      do j=1,3
103         al(j)=0.d0
104         fv1(j)=0.d0
105         fv2(j)=0.d0
106         do i=1,3
107            alz(j,i)=0.d0
108         enddo
109      enddo
110!
111      call rs(ndim,ndim,cmm,al,matz,alz,fv1,fv2,ier)
112!
113!           ------normalizing eigenvectors to P vectors----------
114!
115      dd=dsqrt(alz(1,3)**2+alz(2,3)**2+alz(3,3)**2)
116      p3(1)=alz(1,3)/dd
117      p3(2)=alz(2,3)/dd
118      p3(3)=alz(3,3)/dd
119!
120      do l=1,3
121         do k=1,3
122            cmm(k,l)=0.d0
123            do j=1,3
124               cm(j,l,k)=0.d0
125            enddo
126         enddo
127      enddo
128!
129      do l=1,3
130         do j=1,3
131            do i=1,3
132               do k=1,3
133                  cm(l,j,i)=cm(l,j,i)+c(l,k,j,i)*xn(k);
134               enddo
135            enddo
136         enddo
137      enddo
138!
139      do  l=1,3
140         do  j=1,3
141            do  i=1,3
142               cmm(l,j)=cmm(l,j)+cm(l,j,i)*p3(i);
143            enddo
144         enddo
145      enddo
146!
147      do j=1,3
148         do l=1,3
149            v(j)=v(j)+cmm(l,j)*p3(l)
150         enddo
151      enddo
152!
153      dd=dsqrt(v(1)**2+v(2)**2+v(3)**2)
154      speed=dd/dsqrt(rho*al(3))
155c      speed=dsqrt(dd/rho)
156!
157!     inverse wave speed
158!
159      a=1.d0/speed
160!
161      return
162      end
163
164