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 pk_y0_yg(p2p1,beta,kappa,y0,yg)
20!
21      implicit none
22!
23      real*8 p2p1,beta,kappa,y0,yg,pcrit
24!
25!     adiabatic expansion factor y0 measured (eq.15-17)
26!
27!     author: Yannick Muller
28!
29      pcrit=(2.d0/(kappa+1.d0))**(kappa/(kappa-1.d0))
30
31      if(p2p1.ge.0.63d0) then
32         y0=1d0-(0.41d0+0.35d0*beta**4.d0)/kappa*(1.d0-p2p1)
33      else
34         y0=1d0-(0.41d0+0.35d0*beta**4.d0)/kappa*(1.d0-0.63d0)
35     &        -(0.3475d0+0.1207d0*beta**2.d0-0.3177d0*beta**4.d0)
36     &        *(0.63d0-p2p1)
37!
38      endif
39!
40!     adiabatic expension factor yg isentropic eq 18
41!
42      if(p2p1.ge.1d0) then
43         yg=1.d0
44!
45      elseif (p2p1.ge.pcrit) then
46         yg=p2p1**(1.d0/kappa)*dsqrt(kappa/(kappa-1.d0)
47     &        *(1.d0-p2p1**((kappa-1.d0)/kappa)))/dsqrt(1.d0-p2p1)
48!
49      else
50!     critical pressure ratio
51         yg=(2.d0/(kappa+1.d0))**(1.d0/(kappa-1.d0))
52     &       *dsqrt(kappa/(kappa+1.d0))/dsqrt(1.d0-p2p1)
53      endif
54!
55      return
56!
57      end
58