solve5.f90

运载火箭用拉瓦尔喷管计算（此程序参数保密

subroutine solve5()
use module1
integer eqi
external f
double precision x,f
real ans(3),init(3)
do ii=1,nz
   if(ii<=10)then
   init(1)=tt1*0.88
   init(2)=2.8
   init(3)=pp1+0.01
   else 
   init(1)=tt1*0.90
   init(2)=2.8
   init(3)=pp1*0.9
   end if
   do eqi=1,3
      x=init(eqi)
      eps=1.0E-5
      L=200
      call DATKN(X,EPS,F,L,eqi)
      ans(eqi)=x
   end do
   Te(ii)=ans(1)*T0
   ma(ii)=ans(2)
   Pe(ii)=ans(3)*P0   
end do
end

function F(X,EQI)
use module1
double precision X,F
integer eqi
  select case(eqi)
     case(1)
       F=x+x**(1/(gama-1))*sqrt((1-x))*sqrt((gama+1)/(gama-1))-At/A(ii)*(2/(gama+1))**(1/(gama-1))
     case(2)
       F=x+1/x*(2/(gama+1)*(1+(gama-1)*x*x/2))**((gama+1)/(2*gama-2))
     case(3)
       F=x+x**(1/gama)*sqrt(1-x**((gama-1)/gama))*sqrt((gama+1)/(gama-1))*((gama+1)/2)**(1/(gama-1))-At/A(ii)
     case default
       write(*,*)"error!"
  end select
end function 




subroutine DATKN(X,EPS,F,L,eqi)
double precision X,F,U,V
integer eqi
    M=L
    L=0
10	L=L+1
    U=F(X,eqi)
    V=F(U,eqi)
	if (abs(U-V).LT.EPS*1.0D0)then
	  X=V
	else
	  X=V-(V-U)*(V-U)/(V-2.0*U+X)
	  if (L.NE.M) goto 10
	end if
	return
end subroutine