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

#define SETN       201
#define FREQ_S     4.0e9
#define FREQ_E     6.0e9
#define AA         4.740  /*waveguide's size (cm) */
#define OFFSET     1.875  /*offset length of short circuit 2 (cm) */   
#define PI 3.141592654

typedef struct{ double Re,Im;}complex;

complex com(rc,ic)
double rc,ic;
{
 complex oolong;

 oolong.Re=rc;
 oolong.Im=ic;

 return oolong;
}

double cabs(cx)
complex cx;
{
 double ccx;

 ccx=cx.Re*cx.Re+cx.Im*cx.Im;
 ccx=sqrt(ccx);
 
 return ccx;
}

complex ac(cx,cy)
complex cx,cy;
{
 complex oolong;
 
 oolong.Re=cx.Re+cy.Re;
 oolong.Im=cx.Im+cy.Im;

 return oolong;
}

complex sc(cx,cy)
complex cx,cy;
{
 complex oolong;
 
 oolong.Re=cx.Re-cy.Re;
 oolong.Im=cx.Im-cy.Im;

 return oolong;
}

complex mc(cx,cy)
complex cx,cy;
{
 complex oolong;
 
 oolong.Re=cx.Re*cy.Re-cx.Im*cy.Im;
 oolong.Im=cy.Re*cx.Im+cy.Im*cx.Re;

 return oolong;
}

complex dc(cx,cy)
complex cx,cy;
{
 complex oolong;

 oolong.Re=(cx.Re*cy.Re+cx.Im*cy.Im)/(cy.Re*cy.Re+cy.Im*cy.Im);
 oolong.Im=(cy.Re*cx.Im-cy.Im*cx.Re)/(cy.Re*cy.Re+cy.Im*cy.Im);

 return oolong;
}

double phase(cx)
complex cx;
{
 double pp;
 
 if(cx.Re==0)
  { 
   if(cx.Im>0)
    pp=PI/2.0;
   else if(cx.Im<0)
    pp=PI/(-2.0);
   else
    pp=0.0;
  }
 else
  {
   pp=cx.Im/cx.Re;
   pp=atan(pp);
   if(pp>0)
    {
     if(cx.Re>0)
      pp=pp;
     else
      pp=pp-PI;
    }
   else
    {
     if(cx.Re>0)
      pp=pp;
     else
      pp=pp+PI;
    }
  }
 return pp;
}

/***********************************************************
*******************************************************/
complex short1[SETN],short2[SETN],short22[SETN],ss[8][SETN];
complex edf[SETN],erff[SETN],esf[SETN],s11m[SETN],s11a[SETN];
  
get_cal_data()
{
 FILE *oolong[7];
 char filename[7][15];
 int i,j,ln;
 double ooo;

/* Load the short circuit cal data */
 printf("Input file name for short circuit......\n");
 scanf("%s",filename[0]);
 oolong[0]=fopen(filename[0],"r");
 for(i=0;i<SETN;i++)
  fscanf(oolong[0],"%le %le",&short1[i].Re,&short1[i].Im);
 fclose(oolong[0]);

/* Load the second short circuit cal data */
 printf("Input file name for short circuit 2......\n");
 scanf("%s",filename[1]);
 oolong[1]=fopen(filename[1],"r");
 for(i=0;i<SETN;i++)
  fscanf(oolong[1],"%le %le",&short2[i].Re,&short2[i].Im);
 fclose(oolong[1]);

/* Load the cal data for the load */
/* printf("How many load cal data set of different positions?\n");
 scanf("%d",&ln);*/
 for(i=2;i<3;i++)
  {
   printf("Input file name for the load %d......\n",i-1);
   scanf("%s",filename[i]);
   oolong[i]=fopen(filename[i],"r");
   for(j=0;j<SETN;j++)
    fscanf(oolong[i],
     "%le %le",&ss[i-2][j].Re,&ss[i-2][j].Im);
   fclose(oolong[i]);
  }
}

get_Edf()
{
 int s1,s2,s3,s4,i;
 double aa[SETN],bb[SETN],cc[SETN],dd[SETN];
 
 for(i=0;i<SETN;i++)
  {
/*   if(i<=101)
    {
     s1=0;s2=3;s3=2;s4=4; 
    }
   else
    {
     s1=0;s2=2;s3=1;s4=3;
    }    
   aa[i]=(ss[s2][i].Re-ss[s1][i].Re)/(ss[s1][i].Im-ss[s2][i].Im);
   cc[i]=(ss[s3][i].Re-ss[s4][i].Re)/(ss[s4][i].Im-ss[s3][i].Im);
   bb[i]=0.5*(ss[s1][i].Im+ss[s2][i].Im)
         -0.5*aa[i]*(ss[s1][i].Re+ss[s2][i].Re);
   dd[i]=0.5*(ss[s4][i].Im+ss[s3][i].Im)
         -0.5*cc[i]*(ss[s4][i].Re+ss[s3][i].Re);

   edf[i].Re=(dd[i]-bb[i])/(aa[i]-cc[i]);
   edf[i].Im=(aa[i]*dd[i]-bb[i]*cc[i])/(aa[i]-cc[i]);
*/
   edf[i].Re=ss[0][i].Re;
   edf[i].Im=ss[0][i].Im;
  }
}

get_Erf_Esf()
{
 complex aaa[SETN],bbb[SETN],ccc[SETN],ddd[SETN];
 int i;
 double phase,ww,www;
 
 ww=2.0*PI*FREQ_S;
 www=2.0*PI*(FREQ_E-FREQ_S)/((double)(SETN-1));

 for(i=0;i<SETN;i++)
 {
  phase=-2.0*sqrt(ww*ww/9.0e16-PI*PI/(AA*AA*1.0e-4))*OFFSET*1.0e-2;
  short22[i]=com(-1.0*cos(phase),-1.0*sin(phase));
  printf("%lf %lf %le\n",short22[i].Re,short22[i].Im,ww/(2*PI));
  ww=ww+www;
 }

 for(i=0;i<SETN;i++)
 {
  aaa[i]=sc(short1[i],edf[i]);
  ccc[i]=sc(short2[i],edf[i]);
  bbb[i]=dc(aaa[i],com(-1.0,0.0));
  ddd[i]=dc(ccc[i],short22[i]);

  erff[i]=dc(sc(mc(aaa[i],ddd[i]),mc(bbb[i],ccc[i])),
          sc(aaa[i],ccc[i]));
  esf[i]=dc(sc(bbb[i],ddd[i]),sc(aaa[i],ccc[i]));
 }
}

get_cal()
{
 FILE *oolong;
 char filename[15];
 int i;
 
 printf("Input the filename of data you want to calibrate...\n");
 scanf("%s",filename);

 oolong=fopen(filename,"r");
 for(i=0;i<SETN;i++)
  fscanf(oolong,"%le %le",&s11m[i].Re,&s11m[i].Im);
}

make_cal()
{
 int i;

 for(i=0;i<SETN;i++)
  {
   s11a[i]=dc(sc(s11m[i],edf[i])
           ,ac(mc(esf[i],sc(s11m[i],edf[i])),erff[i]));
  }
}
 
save_cal()
{
 FILE *oolong;
 char filename[15];
 int i;
 double ooop;

 printf("Input the filename of the PHASE of calibrated data...\n");
 scanf("%s",filename);

 oolong=fopen(filename,"w");
 for(i=0;i<SETN;i++)
  { 
   ooop=phase(s11a[i]);
   ooop=ooop*180.0/PI;
   fprintf(oolong,"%le\n",ooop);
  }
 fclose(oolong);

 printf("Input filename of the MAG of calibrated data...\n");
 scanf("%s",filename);

 oolong=fopen(filename,"w");
 for(i=0;i<SETN;i++)
  {
   ooop=cabs(s11a[i]);
   ooop=20.0*log10(ooop);
   fprintf(oolong,"%le\n",ooop);
  }
 fclose(oolong);
}

main()
{
 complex d,d1,d2;

 get_cal_data();
 get_Edf();
 get_Erf_Esf();
 get_cal();
 make_cal();
 save_cal();
}