pt2_lim_calc.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	pt2_lim_calc (pt1, a2, a1, kappa, zeta, pt2_lim)

Function/Subroutine Documentation

pt2_lim_calc()

subroutine pt2_lim_calc	(	real*8	pt1,
		real*8	a2,
		real*8	a1,
		real*8	kappa,
		real*8	zeta,
		real*8	pt2_lim
	)

!
      implicit none
!
      integer i 
!
      real*8 pt1,a2,a1,kappa,pt2_lim,x,zeta,f,df,expon1,
     &     expon2,expon3,cte,a2a1,kp1,km1,delta_x,fact1,fact2,term
!
      x=0.999
!
!     x belongs to interval [0;1]
!
!     modified 25.11.2007 
!     since Pt1/Pt2=(1+0.5(kappa)-M)**(zeta*kappa)/(kappa-1)
!     and for zeta1 elements type M_crit=M1=1
!     and for zeta2 elements type M_crit=M2 =1
!     it is not necessary to iteratively solve the flow equation.
!     Instead the previous equation is solved to find pt2_crit
      if(zeta.ge.0d0) then
         kp1=kappa+1d0
         km1=kappa-1d0
         a2a1=a2/a1
         expon1=-0.5d0*kp1/(zeta*kappa)
         expon2=-0.5d0*kp1/km1
         cte=a2a1*(0.5*kp1)**expon2
         expon3=-km1/(zeta*kappa)
         i=0
!
!        
         do
            i=i+1
!     
            f=x**(-1d0)-cte*x**(expon1)
     &           *(2d0/km1*(x**expon3-1.d0))**-0.5d0
!     
            df=-1.d0/x**2-cte*(x**expon1
     &           *(2d0/km1*(x**expon3-1.d0))**-0.5d0)
     &           *(expon1/x-1d0/km1*expon3*x**(expon3-1d0)
     &           *(2d0/km1*(x**expon3-1.d0))**(-1.d0))
            
            delta_x=-f/df
!     
            if(( dabs(delta_x/x).le.1.e-8)
     &           .or.(dabs(delta_x/1d0).le.1.e-10)) then
!
               pt2_lim=pt1*x
!
               exit
            endif
            if(i.gt.25)then
                pt2_lim=pt1/(1+0.5*km1)**(zeta*kappa/km1)
                exit
             endif
!     
            x=delta_x+x
!     
         enddo
!
      else
!
         do 
            kp1=kappa+1d0
            km1=kappa-1d0
            a2a1=a2/a1
            expon1=kp1/(zeta*kappa)
            expon2=km1/(zeta*kappa)
            expon3=kp1/km1
            cte=a2a1**2*(0.5*kp1)**-expon3*(2/km1)**-1
            fact1=x**-expon1
            fact2=x**-expon2
            term=fact2-1
!     
            f=x**-2-cte*fact1*term**-1
!     
            df=-2*x**-3-cte*(x**(-expon1-1)*term**-1)
     &           *(-expon1+expon2*(x**-expon2)*fact2*term**-1)
!     
            delta_x=-f/df
!     
            if(( dabs(delta_x/x).le.1.e-8)
     &           .or.(dabs(delta_x/1d0).le.1.e-10)) then
               pt2_lim=pt1*x
               exit
            endif
!     
            x=delta_x+x
!     
         enddo
!     
      endif
      
      return