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 calceigenvalues(c,al)
20!
21!     calculates the eigenvalues al of the symmetric 3x3 matrix c
22!     the eigenvalues are sorted in increasing order
23!
24      implicit none
25!
26      integer three,kflag,i
27!
28      real*8 c(3,3),al(3),v1,v2,v3,bb,cc,cm,cn,tt,pi
29!
30      data kflag /1/
31      data three /3/
32!
33!     calculation of the eigenvalues of c
34!     Simo & Hughes Computational Inelasticity p 244
35!
36      pi=4.d0*datan(1.d0)
37!
38      v1=c(1,1)+c(2,2)+c(3,3)
39      v2=c(2,2)*c(3,3)+c(1,1)*c(3,3)+c(1,1)*c(2,2)-
40     &     (c(2,3)*c(2,3)+c(1,3)*c(1,3)+c(1,2)*c(1,2))
41      v3=c(1,1)*(c(2,2)*c(3,3)-c(2,3)*c(2,3))
42     &     -c(1,2)*(c(1,2)*c(3,3)-c(1,3)*c(2,3))
43     &     +c(1,3)*(c(1,2)*c(2,3)-c(1,3)*c(2,2))
44!
45      bb=v2-v1*v1/3.d0
46      cc=-2.d0*v1**3/27.d0+v1*v2/3.d0-v3
47      if(dabs(bb).le.1.d-10) then
48         if(dabs(cc).gt.1.d-10) then
49            al(1)=-cc**(1.d0/3.d0)
50         else
51            al(1)=0.d0
52         endif
53         al(2)=al(1)
54         al(3)=al(1)
55      else
56         cm=2.d0*dsqrt(-bb/3.d0)
57         cn=3.d0*cc/(cm*bb)
58         if(dabs(cn).gt.1.d0) then
59            if(cn.gt.1.d0) then
60               cn=1.d0
61            else
62               cn=-1.d0
63            endif
64         endif
65         tt=datan2(dsqrt(1.d0-cn*cn),cn)/3.d0
66         al(1)=cm*dcos(tt)
67         al(2)=cm*dcos(tt+2.d0*pi/3.d0)
68         al(3)=cm*dcos(tt+4.d0*pi/3.d0)
69      endif
70      do i=1,3
71         al(i)=al(i)+v1/3.d0
72      enddo
73!
74!     sorting
75!
76      call insertsortd(al,three)
77c      call dsort(al,idummy,three,kflag)
78!
79      return
80      end
81