#include "math.h"
#include "stdio.h"
double fd(R,z,d)                                   
double R,z,d;
{
  double result;
  result=6.39*1e8/(350.0*24*3600)-0.03848*sqrt((pow(6.3*1e6,2)
	  -pow(2.0*1e6,2))*pow(d,5.0)/(R*z*0.5624*300.15*208000.0));
  return (result);
}
  double fr(k,d,Re,Rr)
  double k,d,Re,Rr;
{
	double result;
	result=1/sqrt(Rr)+2*log10(k/(3.7*d)+2.51/(Re*sqrt(Rr)));
	return (result);
}
void main()
{ 
    double u,x,y,cu,z,Pq,Pz;
    double Q,Pi,T,Ro,Pr,Tr,Pc,Tc,Ppj;
    double Re,R,R1,R2,R3,d,d1,d2,d3;
    double del,c,k;
      c=0.03848;  del=0.5624;  T=300.15;  k=0.05*1e-3; 
      Q=6.39*1e8/(350*24*3600);  Pq=6.3*1e6;  Pz=2*1e6;
      Pi=3.1415926;   Ro=0.7272;
      Pc=4.625*1e6;   Tc=191.3;
      x=2.57+0.2781*del+1063.6/T;
      y=1.11+0.04*x;
      cu=2.415*(7.77+0.1844*del)*pow(T,1.5)/(122.4+377.58*del+1.8*T)*1e-4;
      u=cu*exp(x*pow((Ro/1000),y));									        //計(jì)算動(dòng)力粘度
      Ppj=2.0/3.0*(Pq+Pz*Pz/(Pq+Pz));								        //計(jì)算平均壓力
      Pr=Ppj/Pc;													        //計(jì)算對(duì)比壓力
      Tr=T/Tc;                                                              //計(jì)算對(duì)比溫度
      z=Pr*(0.1391*Tr-0.2988)+0.0007*Tr+0.9969;                             //計(jì)算壓縮系數(shù)
    printf ("對(duì)比壓力為Pr=%e,對(duì)比溫度為Tr=%e\n\n\n",Pr,Tr);
    printf ("動(dòng)力粘度為u=%e毫帕.秒\n\n\n",u);
    printf ("氣體平均壓縮系數(shù)為Z=%e\n\n\n",z);
      d1=0.001;d2=1.5;
      while(1)
{
	   d3=(d1+d2)/2.0;
       Re=1.534*Q*del*1000/(d3*u);
       R1=0.001;
       R2=1.0;
      while(1)
	  {
	   R3=(R1+R2)/2.0;
	   if(fr(k,d3,Re,R3)>1e-6) R1=R3;
	   if(fr(k,d3,Re,R3)<-1e-6) R2=R3;
	   if(fabs(fr(k,d3,Re,R3))<1e-6) break;
	  }
       R=R3;
       if (fd(R,z,d3)>1e-6) d1=d3;
       if (fd(R,z,d3)<-1e-6) d2=d3;
       if (fabs(d3-d1)<1e-6&&fabs(d3-d2)<1e-6) break;
}
       d=d3;
     printf("摩阻系數(shù)為R=%e,管徑為d=%em\n\n\n",R,d);
     printf("雷諾數(shù)為Re=%e\n\n\n",Re);
}