#include <stdio.h>
#include <math.h>

#define epsilon 0.000001f
#define shag 0.05f

double sqr(double p)
{
  return p*p;
}

double Y(double y, double a[])
{
  return (double)(a[0]*pow(y,6)+a[1]*pow(y,5)+a[2]*pow(y,4)+a[3]*pow(y,3)+a[4]*pow(y,2)+a[5]*y+a[6]);
}

double mod(double x)
{
  if (x>0.0f) return x; else return (-x);
}

double solve_range(double a[])
{
  double A,B,max,min;
  int i;
  A=0.0f;B=0.0f;
  for(i=0;i<7;i++)
    if (mod(a[i])>A) A=mod(a[i]);
  for(i=0;i<6;i++)
    if (mod(a[i])>B) B=mod(a[i]);
  max=1+A/mod(a[0]);
  min=mod(a[6])/(mod(a[6])+B);
  printf("Solve range:nmax %lf nmin %lf nn",max,min);
  return max;
}

double halfsolve(double a[],double aa,double bb)
{
  double cc;
  cc=(aa+bb)/2;
  if (mod(bb-aa)>epsilon)
  {
  if (Y(aa,a)*Y(cc,a)<0.0f)
      return halfsolve(a,aa,cc); else return halfsolve(a,cc,bb);
  } else return cc;  
}

int solve(double a[],double range,double sol[])
{
  int i;
  double c,b,solv;
  i=0;
  c=0.0f;
  while (b<range)
    {
     b=c;
     c+=shag;
     if (Y(c,a)*Y(b,a)<0.0f) 
       {
 printf("local %lf %lfn",b,c);
 solv=halfsolve(a,b,c);
 //printf("solve %lfn",solv);
 sol[i]=solv;
 i++;
       }
    }
    return (i-1);
}

int main()
{
double l0,p0,P0,U0;
double l3,p3,U3,P3;

double a[7],Y[7],j0,j3,X;
double C0,C3,e0,e3;
int i,solvz;
double P1,p1,p2,D0,D3,U1;

l0=5/3;
p0=2.71*pow(10,-3);
P0=5*pow(10,5);
U0=0.0f;

l3=7/5;
p3=2.219*pow(10,-4);
U3=1.587*pow(10,5);
P3=pow(10,6);

j0=(l0+1)/(l0-1);
j3=(l3+1)/(l3-1);
X=P3/P0;

C0=sqrt(l0*P0/p0);
C3=sqrt(l3*P3/p3);

printf("C0 %lf nC3 %lf nn",C0,C3);

e0=2*pow(C0,2)/(l0*(l0-1)*sqr(U3-U0));
e3=2*pow(C3,2)/(l3*(l3-1)*sqr(U3-U0));

printf("e0 %lf ne3 %lf nn",e0,e3); 

a[0]=sqr(j0*e3-j3*X*e0);
a[1]=2*((j0*e3-j3*X*e0)*(e3*(1-2*j0*X)-e0*X*(X-2*j3))-j0*j3*X*(j0*e3+j3*X*e0));
a[2]=sqr(e3)*(6*sqr(j0)*sqr(X)-8*j0*X+1)-2*e0*e3*X*(j0*j3*(X*X+4*X+1)-2*(X+1)*(j3+j0*X)+X)+sqr(e0)*sqr(X)*(6*sqr(j3)-8*j3*X+X*X)+sqr(j0*j3*X)-2*j0*X*e3*(j0*X-2*j0*j3*X+2*j3)-2*j3*X*X*e0*(j3+2*j0*X-2*j0*j3);
a[3]=-2*X*(2*e3*e3*(sqr(j0*X)-3*j0*X+1)+e0*e3*((j3+j0*X)*(X*X+4*X+1)-2*j0*j3*X*(X+1)-2*X*(X+1))+2*e0*e0*X*(X*X-3*j3*X+j3*j3)-j0*j3*X*(j0*X+j3)+e3*(j0*j0*j3*X*X-2*X*(2*j0*j3+j0*j0*X)+(2*j0*X+j3))+e0*X*(j0*j3*j3-2*j3*(j3+2*j0*X)+2*j3*X+j0*X*X));
a[4]=X*X*(e3*e3*(sqr(j0*X)-8*j0*X+6)-2*e0*e3*(j0*j3*X-2*(X+1)*(j3+j0*X)+X*X+4*X+1)+e0*e0*(j3*j3-8*j3*X+6*X*X)+(j3*j3+4*j0*j3*X+sqr(j0*X))-2*e3*((j0*j0*X+2*j0*j3)*X-2*(2*j0*X+j3)+1)-2*e0*(j3*(2*j0*X+j3)-2*X*(2*j3+j0*X)+X*X));
a[5]=2*X*X*X*(e3*e3*(j0*X-2)-e0*e3*(j0*X-2+j3-2*X)+e0*e0*(j3-2*X)+(j3+j0*X)-e3*(2*j0*X+j3-2)-e0*(2*j3+j0*X-2*X));
a[6]=pow(X,4)*(sqr(e3-e0)+1-2*(e3+e0));

for (i=0;i<7;i++)
   printf("a%i %f n",i,a[i]);
printf("n");

solvz=solve(a,solve_range(a),Y);

for (i=0;i<solvz+1;i++)
{
   printf("Y%i %f n Solv D for itn",i,Y[i]);
   P1=Y[i]*P0;
   p1=p0*((l0-1)+(l0+1)*P1/P0)/((l0+1)+(l0-1)*P1/P0);
   U1=U0-C0*sqrt(2)*(Y[i]-1)/(sqrt(l0*(l0-1))*sqrt(1+j0*Y[i]));
   D0=(p0*U0-p1*U1)/(p0-p1);
   p2=p3*((l3-1)+(l3+1)*P1/P3)/((l3+1)+(l3-1)*P1/P3);
   D3=(p2*U1-p3*U3)/(p2-p3);
   printf(" P1 %lfn p1 %lfn U1,2 %lfn D0 %lfn",P1,p1,U1,D0);
   printf(" p2 %lfn D3 %lfn",p2,D3);
}
printf("n");

return 0;
}