#include"stdlib.h"
#include"math.h"
#include"conio.h"
#include"stdio.h"
#include<string.h>
#define N 100
#define M 15
#define Ur 9.524
#define W  20
#define PAI 3.14

void chazhi_1(x,y,g,s,dfw)
float x[M],y[M],g[N],s[N],dfw[N];
{
 float h[M],m[M],a[M],b[M],c[M],d[M],p[M],q[M];
 int i,j,k;
 for(i=0;i<M;i++)
 b[i]=2;
 for(j=1;j<M;j++)
 h[j]=x[j]-x[j-1];
 m[0]=0;m[M-1]=0;d[0]=2*m[0];d[M-1]=2*m[M-1];c[0]=0;a[M-1]=0;c[M-1]=0;
 for(j=1;j<M-1;j++)
 {
  c[j]=h[j+1]/(h[j]+h[j+1]);
  a[j]=1-c[j];
  d[j]=6/(h[j]+h[j+1])*((y[j+1]-y[j])/h[j+1]-(y[j]-y[j-1])/h[j]);
 }
  q[0]=c[0]/b[0];p[0]=d[0]/b[0]; /*q=c  p=d*/
  for(i=1;i<M;i++)
 {
  q[i]=c[i]/(b[i]-a[i]*q[i-1]);
  p[i]=(d[i]-a[i]*p[i-1])/(b[i]-a[i]*q[i-1]);
 }
  m[M-1]=p[M-1];
  for(i=M-2;i>-1;i--)
  m[i]=p[i]-q[i]*m[i+1];
  for(j=0;j<N;j++)
  for(i=1;i<M;i++)
  if(g[j]<=x[i])
 {
  s[j]=m[i-1]*pow((x[i]-g[j]),3)/h[i]/6+m[i]*pow((g[j]-x[i-1]),3)/h[i]/6+(y[i-1]-m[i-1]*h[i]*h[i]/6)*(x[i]-g[j])/h[i]+(y[i]-m[i]*h[i]*h[i]/6)*(g[j]-x[i-1])/h[i];
  dfw[j]=-m[i-1]*pow((x[i]-g[j]),2)/2/h[i]+m[i]*pow((g[j]-x[i-1]),2)/2/h[i]+(y[i]-y[i-1])/h[i]-h[i]/6*(m[i]-m[i-1]);
  break;
 }
}

void chazhi_2(x,y,g,s)
float x[M],y[M],g[N],s[N];
{
 float h[M],m[M],a[M],b[M],c[M],d[M],p[M],q[M];
 int i,j,k;
 for(j=1;j<M;j++)
 h[j]=x[j]-x[j-1];
 m[0]=0;m[M-1]=0;d[0]=2*m[0];d[M-1]=2*m[M-1];c[0]=0;
 for(j=1;j<M-1;j++)
 {
  c[j]=h[j+1]/(h[j]+h[j+1]);
  a[j]=1-c[j];
  d[j]=6/(h[j]+h[j+1])*((y[j+1]-y[j])/h[j+1]-(y[j]-y[j-1])/h[j]);
 }
 for(i=0;i<M;i++)
 b[i]=2;
 q[0]=c[0]/b[0];p[0]=d[0]/b[0];
 for(i=1;i<M;i++)
 {
  q[i]=c[i]/(b[i]-a[i]*q[i-1]);
  p[i]=(d[i]-a[i]*p[i-1])/(b[i]-a[i]*q[i-1]);
 }
 m[M-1]=p[M-1];
 for(i=M-2;i>-1;i--)
 m[i]=p[i]-q[i]*m[i+1];
 for(j=0;j<N;j++)
 for(i=1;i<M;i++)
 if(g[j]<=x[i])
 {
  s[j]=m[i-1]*pow((x[i]-g[j]),3)/h[i]/6+m[i]*pow((g[j]-x[i-1]),3)/h[i]/6+(y[i-1]-m[i-1]*h[i]*h[i]/6)*(x[i]-g[j])/h[i]+(y[i]-m[i]*h[i]*h[i]/6)*(g[j]-x[i-1])/h[i];
  break;
 }
}
float chazhi(float x1,float z1,float x2,float z2,float x)
{
 float z;
 z=z1+(x-x1)*(z2-z1)/(x2-x1);
 return(z);
}

main()
{
 float chazhi_();
 void  chazhi_1();
 void  chazhi_2();
 float mid,swf,dsw;
 float sw[N],kro[N],krw[N],fw[N],dfw[N];
 float sw0[M],fw0[M],krw0[M],kro0[M];
 float re,qo11,t2i_,t2b,tjb,t2b_,t2off;
 float ta[N],ta_[N],qoo[N],qto[N],qcw[N],caishoulv[N];
 float kroswc,aa,bb,kk,pp,qoo1[N],qoo2[N],hw[N],qo[N],qw[N];
 float e,t1b_,dt,t_[N],q_[N],tb_[N],qlt_[N],qot_[N],qwt_[N],fwe[N],qo2[N],swe[N],sw_ave[N];
 float pw,pwf1,pwf2,vp1,vp2,ss[N],z[N];
 float B,H,K,MD,U,L1,L2,Nw,No1,No2,do1,do2,re1,re2,ri1,ri2;
 float t1b,qi1,qi2,Qi1,Qit,Qi2 ,swc,v1,v2,t1off,q2i,F,dfwswf,t1off_,nf[N],Q1;
 float riw,Bo,Rw,Ro;
 FILE *fp;
 int i,j,y1,y2,a,y3,y4,k,b;
 float  k1,k2 ;
 clrscr();
 if((fp=fopen("d:\\dljichu.txt","r"))==NULL)
 {
  printf("dljichu cannot open file\n");
  exit(1);
 }
 fscanf(fp,"%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f",&B,&H,&F,&K,&MD,&U,&Bo,&L1,&L2,&Nw,&No1,&No2,&Rw,&Ro,&pw,&pwf1,&pwf2,&fw);
/*B5000 H12 F0.22 K0.28 MD855 U0.01 Bo1.15 L1950 L2450 Nw20 No115 No215 Rw0.12 Ro0.14 Pw26 Pwf116 Pwf26 fw0.86*/
 fclose(fp);
 if((fp=fopen("d:\\dlswkrow.txt","r"))==NULL)
 {
  printf("dlswkrow.dat cannot open file\n");
  exit(1);
 }
 for(i=0;i<M;i++)
 fscanf(fp,"%f%f%f",&sw0[i],&kro0[i],&krw0[i]);
 fclose(fp);
 dsw=(sw0[M-1]-sw0[0])/(N-1);
 kroswc=kro0[0];
 sw[0]=sw0[0];
 swc=sw[0];
 for(i=1;i<N;i++)
 sw[i]=sw[i-1]+dsw;
 fw0[0]=0.000;
 for(i=1;i<M;i++)
 fw0[i]=1.0/(1.0+kro0[i]/(krw0[i]*Ur));
 chazhi_1(sw0,fw0,sw,fw,dfw);
 chazhi_2(sw0,kro0,sw,kro);
 chazhi_2(sw0,krw0,sw,krw);
 mid=(fw[1]-fw[0])/(sw[1]-sw[0]);
 k=i;
 for(i=1;i<N;i++)
 {
  mid=(fw[i]-fw[0])/(sw[i]-sw[0]);
  if(dfw[i]>mid)
  {
   mid=dfw[i];k=i;
  }
 }
 printf("k=%d\n",k);
 dfwswf=dfw[k];
 swf=sw[k];
 if((fp=fopen("d:\\dlcanshu.txt","w"))==NULL)
 {
  printf("dlcanshu cannot open file\n");
  exit(1);
 }
 fprintf(fp,"     sw         krw          kro        fw         dfw\n");
 for(i=0;i<N;i++)
 fprintf(fp,"%11.5f%11.5f%11.5f%11.5f%11.5f\n",sw[i],krw[i],kro[i],fw[i],dfw[i]);
 fprintf(fp,"k=%2d,swf=%0.5f\n",k,swf);
 fclose(fp);
 riw=0.5*U*log(0.5*B/Nw/(PAI*Rw))/(Nw*kroswc*K*pow(10,-12)*PAI*H);
 re1=U*L1/(B*kroswc*H*K*pow(10,-12));
 ri1=0.5*U*log(0.5*B/No1/PAI/Ro)/No1/(K*pow(10,-12))/kroswc/PAI/H;
 ri2=0.5*U*log(0.5*B/No2/PAI/Ro)/No2/(K*pow(10,-12))/kroswc/PAI/H;
 re2=U*L2/(B*K*pow(10,-12)*kroswc*H);
 k1=(pwf1-pwf2)*(riw+re1)/(re2+ri2);
 k2=(riw+re1)*ri1/(re2+ri2);
 Qi1=(pw-pwf1-k1)/(riw+re1+ri1+k2)*60*60*365*24;
 Q1=(pw-pwf1-k1)/(riw+re1+ri1+k2);
 qi1=Qi1/No1;
 Qi2=(pwf1-pwf2+Q1*ri1)/(re2+ri2)*60*60*365*24;
 printf("Qi1=%f\nQi2=%f\n",Qi1,Qi2);
 qi2=Qi2/No2;
 Qit=Qi1+Qi2;
 vp1=B*H*F*L1;
 vp2=B*H*F*L2;
 v1=vp1*Qi1/Qit;
 v2=vp2+vp1*Qi2/Qit;
 printf("v1=%f\nv2=%f\nvp1=%f\nvp2=%f\n",v1,v2,vp1,vp2);
 printf("dfwswf=%f\nkroswc=%f\n",dfwswf,kroswc);
 z[49]=0.0;
 for(i=N-2;i>=k;i--)
 {
  kk=dfw[i]-dfw[i+1];
  aa=kro[i]+Ur*krw[i]+Ur*krw[i+1]+kro[i+1];
  z[i]=2*kk/aa+z[i+1];
 }
 printf("z[k]=%f\n",z[k]);
 e=kroswc*z[k]-dfwswf;
 t1b_=0.5*(2*dfwswf+e)/(dfwswf*dfwswf);
 printf("e=%f,t1b_=%f\n",e,t1b_);
 t1b=t1b_*v1/Qi1;
 dt=t1b_/(W-1);
 t_[0]=0.0;
 q_[0]=0.0;
 fwe[0]=0.0;
 caishoulv[0]=0.0;
 for(i=1;i<W;i++)
 {
  t_[i]=t_[i-1]+dt;
  q_[i]=1/sqrt(1+2*e*t_[i] );
  fwe[i]=0.0;
  caishoulv[i]=(sqrt(1+2*e*t_[i])-1)/e/(1-swc);
 }
 if((fp=fopen("d:\\dlbt1.txt","w"))==NULL)
 {
  printf("dlbt1 cannot open file\n");
  exit(1);
 }
 fprintf(fp,"jianshuiqianwuyincishijian wuyincichanliang hanshuilv caishoulv\n");
 for(j=0;j<W;j++)
 fprintf(fp,"%10.2f      %10.2f        %0.2f      %0.2f\n",t_[j],q_[j],fwe[j],caishoulv[j]);
 fclose(fp);
 ss[k]=0;

 for(i=k+1;i<N-1;i++)
 {
  kk=dfw[i]-dfw[i-1];
  pp=pow(2/(dfw[i]+dfw[i-1]),3);
  bb=z[i]+z[i-1];
  ss[i]=bb/2*pp*kk+ss[i-1];
  tb_[i]=t1b_-kroswc*ss[i];
  fwe[i]=fw[i];swe[i]=sw[i];
  qlt_[i]=dfw[i]/(kroswc*z[i]);
  qot_[i]=qlt_[i]*(1.0-fwe[i]);
  qwt_[i]=qlt_[i]*fwe[i];
  hw[i]=qwt_[i]/qlt_[i];
  sw_ave[i]=(1-fwe[i])/dfw[i]+swe[i];
  caishoulv[i]=(sw_ave[i]-swc)/(1-swc);
 }

 if((fp=fopen("d:\\dlbt2.txt","w"))==NULL)
 {
  printf("dlbt2 cannot pen file\n");
  exit(1);
 }

 fprintf(fp,"jianshuihoudianshu wuyincishijian hanshuilv wuyincichanyeliang wuyincichanyouliang wuyincichanshuiliang hanshuilv caishoulv\n");
 for(i=k+1;i<N;i++)
 fprintf(fp,"%3d     %10.5f    %10.5f    %10.5f    %10.5f   %10.5f    %10.5f   %10.5f\n",i,tb_[i],dfw[i],qlt_[i],qot_[i],qwt_[i],hw[i],caishoulv[i]);
 fclose(fp);

 for(i=k;i<N-1;i++)
 if(fw[i]>0.82)
 {
  t1off_=chazhi(fw[i-1],tb_[i-1],fw[i],tb_[i],0.82);
  break;
 }
 printf("t1off_=%f\n",t1off_);
 t1off=t1off_*v1/Qi1;
 re=U*(L1+L2)/(B*H*K*pow(10,-12));
 qo11=(pw-pwf2)/(riw+re+ri2)*60*60*24*365;
 for(i=k;i<N-1;i++)
 if(fw[i]>0.98)
 {
  tjb=chazhi(fw[i-1],tb_[i-1],fw[i],tb_[i],0.98);
  break;
 }
 printf("tjb=%f\n",tjb);

 t2i_=Qit*t1off/(vp1+vp2);
 t2b=(t1b_-t2i_)*(vp1+vp2)/qo11+t1off;
 t2off=(tjb-t2i_)*(vp1+vp2)/qo11+t1off;
 printf("t2off=%f\nt2b=%f\n",t2off,t2b);
 printf("t1b=%10.3f\nt1off=%10.3f\n",t1b,t1off);
 if((fp=fopen("d:\\dljieguo.txt","w"))==NULL)
 {
  printf("dljieguo cannot open file\n");
  exit(1) ;
 }
 fprintf(fp,"yierjieduanchushichanliangjijianshuishijian\n");
 fprintf(fp,"Qi1=%f\nQi2=%f\nQit=%f\nt1b=%f\nt1off=%f\ne=%f\n",Qi1,Qi2,Qit,t1b,t1off,e);
 fprintf(fp,"qo11=%10.3f\nt2i_=%10.3f\nt2b=%10.3f\nt20ff=%10.3f\n",qo11,t2i_,t2b,t2off);
 fclose(fp);
 for(i=0,j=0;i<=t2off;i++)
 {
  nf[j]=i;j+=1;
  if(i<t1b&&t1b<i+1)
  {
   nf[j]=t1b;y1=j;j=j+1;
  }
  if(i<t1off&&t1off<i+1)
  {
   nf[j]=t1off;y2=j;j=j+1;
  }
  if(i<t2b&&t2b<i+1)
  {
   nf[j]=t2b;y3=j;j=j+1;
  }
  if(i<t2off&&t2off<i+1)
  {
   nf[j]=t2off;y4=j;j=j+1;
  }
 }



 for(i=0;i<=y1;i++)
 ta[i]=nf[i];
 qto[0]=0.0;
 qcw[0]=0.0;
 hw[0]=0.0;
 qoo1[0]=Qi1/sqrt(1+2*e*ta[0]*Qi1/v1);
 qoo2[0]=Qi2/sqrt(1+2*e*ta[0]*Qi2/v2);
 qoo[0]=qoo1[0]+qoo2[0];
 for(j=1;j<=y1;j++)
 {
  qoo1[j]=Qi1/sqrt(1+2*e*ta[j]*Qi1/v1);
  qoo2[j]=Qi2/sqrt(1+2*e*ta[j]*Qi2/v2);
  qoo[j]=qoo1[j]+qoo2[j];
  qcw[j]=0.0;
  hw[j]=0.0;
  qto[j]=qto[j-1]+qoo[j];
 }
 for(i=y1+1;i<=y2;i++)
 {
  ta[i]=nf[i];
  ta_[i]=ta[i]*Qi1/v1;
  for(a=k;a<N;a++)
  if(ta_[i]<=tb_[a])
  {
   qoo1[i]=chazhi(tb_[a-1],qot_[a-1],tb_[a],qot_[a],ta_[i])*Qi1;
   qoo2[i]=chazhi(tb_[a-1],qot_[a-1],tb_[a],qot_[a],ta_[i])*Qi2;
   qoo[i]=chazhi(tb_[a-1],qot_[a-1],tb_[a],qot_[a],ta_[i])*Qit;
   qcw[i]=chazhi(tb_[a-1],qwt_[a-1],tb_[a],qwt_[a],ta_[i])*Qit;
   qto[i]=qto[i-1]+qoo[i];
   hw[i]=qcw[i]/(qcw[i]+qoo[i]) ;
   break;
  }
 }
 for(i=0;i<y2+1;i++)
 caishoulv[i]=qto[i]/(vp1+vp2)/(1-swc);
 for(j=y2+1;j<y3+1;j++)
 {
  ta_[j]=t2i_+qo11*(nf[j]-t1off)/(vp1+vp2);
  qo2[j]=qo11/sqrt(1+2*e*ta_[j]);
  qcw[j]=0.0;
  hw[j]=qcw[j]/(qcw[j]+qo2[j]);
  qto[j]=qto[j-1]+qo2[j];
 }
 for(i=y3+1;i<=y4;i++)
 {
  ta[i]=nf[i];
  ta_[i]=(ta[i]-t1off)*qo11/(vp1+vp2)+t2i_;
  for(b=k+1;b<N;b++)
  if(ta_[i]<=tb_[b])
  {
   qo2[i]=chazhi(tb_[b-1],qot_[b-1],tb_[b],qot_[b],ta_[i])*qo11;
   qcw[i]=chazhi(tb_[b-1],qwt_[b-1],tb_[b],qwt_[b],ta_[i])*qo11;
   hw[i]=qcw[i]/(qcw[i]+qo2[i]);
   qto[i]=qo2[i]+qto[i-1];
   break;
  }
 }
 for(i=y2;i<y4;i++)
 caishoulv[i]=qto[i]/((vp1+vp2)*(1-swc));
 if((fp=fopen("d:\\dlquannian.txt","w"))==NULL)
 {
  printf("dlquannian cannot open file\n");
  exit(1);
 }
 qto[0]=0.0;qcw[0]=0.0;
 for(i=1;i<=y2;i++)
 {
  qo[i]=qto[i]-qto[i-1];qw[i]=qcw[i]-qcw[i-1];
 }
 qw[y2+1]=0.0;qo[y2+1]=qto[y2+1]-qto[y2];
 for(i=y2+2;i<=y4;i++)
 {
  qo[i]=qto[i]-qto[i-1];qw[i]=qcw[i]-qcw[i-1];
 }

 fprintf(fp,"nianfen nianchanyouliang nianchanshuiliang leijichanyouliang leijichanshuiliang caishoulv hanshuilv\n");
 for(i=0;i<=y4;i++)
 fprintf(fp,"%5.2f  %12.3f  %12.3f  %12.3f %10.3f %12.3f %12.3f\n",nf[i],qo[i],qw[i],qto[i],qcw[i],caishoulv[i],hw[i]);
 fclose(fp);
}