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

extern void upup_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
extern void downdown_(int *L1,int *L2, int *L3,int *L4,double *s,double *t,double *u,int *iz,double *har,double *hai);
extern void flagset_(int *iqedx,int *iewx,int *ibornx,int *gfschemex,int *ifggx,double *ncx,double *fcx,double *tlmu2x);
extern void paraget_(double  *mtax, double  *conhcx, double *pix);
extern void printconsts_(int *pmg);


int iqed=0,iew=1,iborn=0,gfscheme=0,ifgg=1;
double nc=1.0,fc=3.0,tlmu2=1e-5;
int buf=0;
double R[4][4];
double divv;
double borndivv=123.456;

struct compl {double im,re;};

struct compl mul(struct compl x, struct compl y)
{
  struct compl c ;
  c.re=x.re*y.re-x.im*y.im;
  c.im=x.re*y.im+x.im*y.re;
  return c;
}

struct compl conju(struct compl x)
{
  struct compl c ;
  c.re=x.re;
  c.im=-x.im;
  return c;
}

double abso(struct compl x)
{
  return sqrt(x.re*x.re+x.im*x.im);
}


void Rcalc(int flav, double sloop,double costhetloop)
{
	int L1,L2,L3,L4,iz,i,j,M3,M4;
	double s,t,u;
	double cosf,betaf,har,hai,MM,sinf,xnorm,xx;

	double mta,conhc,pi;

	double sum,sig[4][4][2], Bench[4][4]={0};
	double ha[2];
      	struct compl amp[2][2][2][2][2]={0}, sigma[4][2][2]={0}, c;
        sigma[0][0][0].re=1;
        sigma[0][1][1].re=1;

        sigma[1][0][1].re=1;
        sigma[1][1][0].re=1;

        sigma[2][0][1].im=-1;
        sigma[2][1][0].im= 1;

        sigma[3][0][0].re=1;
        sigma[3][1][1].re=-1;

	s=sloop;
	cosf=costhetloop;
// --------------------------
	paraget_(&mta,&conhc,&pi);
	betaf = sqrt(1e0-4e0*mta*mta/s);
	t = mta*mta - s/2*(1e0-betaf*cosf);
	u = mta*mta - s/2*(1e0+betaf*cosf);
	for(iz = 0;iz<2;iz++)
	{
		for(L1=1;L1<=2;L1++)
		{
			for(L2=1;L2<=2;L2++)
			{
				for(L3=1;L3<=2;L3++)
				{
					for(L4=1;L4<=2;L4++)
					{
					  if(flav==1)	downdown_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
					  else  	upup_(&L1,&L2,&L3,&L4,&s,&t,&u,&iz,&har,&hai);
                                                amp[L1-1][L2-1][L3-1][L4-1][iz].im=hai;
                                                amp[L1-1][L2-1][L3-1][L4-1][iz].re=har;
					}
				}
			}
		}
	}
	for(i=0;i<4;i++){
	for(j=0;j<4;j++){
	for(iz = 0;iz<2;iz++)
	{
	        sum = 0e0;
		for(L1=1;L1<=2;L1++)
		{
			for(L2=1;L2<=2;L2++)
			{
				for(L3=1;L3<=2;L3++)
				{
					for(L4=1;L4<=2;L4++)
	                                for(M3=1;M3<=2;M3++)
				        for(M4=1;M4<=2;M4++)
					{       
					  c=mul(mul(mul(amp[L1-1][L2-1][L3-1][L4-1][iz],conju(amp[L1-1][L2-1][M3-1][M4-1][iz])),
						    sigma[i][L3-1][M3-1]),(sigma[j][M4-1][L4-1])); 
					  sum=sum+ c.re;
					}
				}
			}
		}
		sig[i][j][iz]=sum;
	}
	R[i][j] = conhc *                            // to pbarn
                  nc/fc*1e0/2.0/s *
	          1e0/4 *                            // spin sum
                  (sig[i][j][1] - sig[i][j][0]) *    // |Amp|^2 - linearized
	          betaf/16/pi;                       // phase_space/dcos{theta}
        if(i==0 && j==0)
	    divv=R[0][0];

	R[i][j]=R[i][j]/divv;  
	        if(i!=0) R[i][j] =-  R[i][j]; // necessary for tau- but further investigation needed

	}
	}

	/*	// there should be rotation like this one, however we do not need that. 
                // Wrong definition of Pauli matrix for antiparticle, missing conjugation?
	for(j=0;j<4;j++)
	  {
	    xx=       R[1][j];
           R[3][j] =-  R[3][j];
           R[1][j] =  -xx;
	  }
	
	*/
	
						
	for(i=0;i<4;i++)      // rotations  around z axis  must be checked
	{
	  xx=       R[i][1];
          R[i][1] =-  R[i][2];
          R[i][2] =  xx;
	}
	
	
	for(j=0;j<4;j++)    // rotations  around z axis  must be checked
        {
	  xx=       R[1][j];
          R[1][j] =  R[2][j];
          R[2][j] =-  xx;
	}
		
				
	/*	 	
	  printf("\n");
          printf("d_sigma/d_cos{theta} = %.9f \n",divv);
	  printf("\n");
	  printf("R(i,j)= \n");
	  printf("\n");
          for(i=0;i<4;i++)
	  {
	    for(j=0;j<4;j++)
	      printf("%.5f ",R[i][j]);
	    printf("\n");
	  }
	 
	  sinf=sqrt(1-cosf*cosf);
          MM=sqrt(4e0*mta*mta/s);
          xnorm=1+cosf*cosf + MM*MM*sinf*sinf;
          Bench[0][0]=(1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
          Bench[1][1]=(-(1- MM*MM)*sinf*sinf)/xnorm;
          Bench[2][2]=( (1+ MM*MM)*sinf*sinf)/xnorm;
          Bench[3][3]=(1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
          Bench[2][3]=(2*MM*sinf*cosf)/xnorm;
          Bench[3][2]=(2*MM*sinf*cosf)/xnorm;
	  printf("\n");
	  printf("Bench(i,j)= \n");
	  printf("\n");

	  for(i=0;i<4;i++)
	  {
	    for(j=0;j<4;j++)
	      printf("%.5f ",Bench[i][j]);
	    printf("\n");
	  }

	  printf("\n");
	  printf("   \n");
	  printf("Bench(i,j)=R(i,j) at low energies, where ew corr and Z contribute little \n");
	  printf("\n");
	*/	  

}

int main(int argc, char **argv){
  double sloop=100;
  double costhetloop=0.0;
  int i,j,k,l;
  int NS1=100, NS2=100, NS3=100, NCOS=21;
  double smin1=log(6*6);
  double smax1=log(17000*17000);

  double smin2=85*85;
  double smax2=110*110;

  double smin3=160*160;
  double smax3=220*220;
  int    flav;
  double step;
  FILE *f=NULL,*d=NULL,*dd=NULL;
  unsigned char c;
  char   ftime[255],path[255];
  time_t utime;
  if (argc > 1) flav=atoi(argv[1]);
  else flav=1;
  if (flav==1) f=fopen("table1-1.txt","w");
  else         f=fopen("table2-2.txt","w");

  flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
  printconsts_(&buf);
  dd=fopen("SancLib_v1_02/lib.txt","r");
  if(!dd) { printf("No info file  on electrowek library SANC (file SancLib_v1_02/lib.txt missing), we better stop \n"); return -1; }
  d=fopen("var.dump","r");
  if(!d) { printf("No initialization variables of SANC (file var.dump missing), we better stop \n"); return -1; }

  fprintf(f,"Dimensions: %i %i %i %i\n",NS1,NS2,NS3,NCOS);
  fprintf(f,"Ranges: %.5f %.5f %.5f %.5f %.5f %.5f \n",smin1,smax1,smin2,smax2,smin3,smax3);
  getcwd(path,255);
  time(&utime);
  strftime(ftime,255,"%d %b %Y, %H:%M:%S, GMT%z",localtime(&utime));
  fprintf(f,"Timestamp: %s\n",ftime);
  fprintf(f,"Path: %s\n\n",path);

  while(1)
  {
    c=0;
    fscanf(dd,"%c",&c);
    if(c==0) break;
    fprintf(f,"%c",c);
    }

  while(1)
  {
    c=0;
    fscanf(d,"%c",&c);
    if(c==0) break;
    fprintf(f,"%c",c);
    }
  fprintf(f,"\nBeginRange1\n");  
  printf("\nBeginRange1\n");  

  step=(smax1-smin1)/(NS1-1);
  for(i=0;i<NS1;i++)
  {
    sloop=exp(smin1+i*step);
    printf("%.5f \n",sloop);
    for(j=0;j<NCOS;j++)
    {
      costhetloop=-1.+1.0/NCOS+j*2./(NCOS);

      iew=0;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);
      borndivv=divv;

      iew=1;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);

      for(k=0;k<4;k++)
	for(l=0;l<4;l++)
	  fprintf(f,"%.5f ",R[k][l]);
      fprintf(f,"%.8f ",divv);
      fprintf(f,"%.8f\n",borndivv);
      // fprintf(f,"%.5f \n",sloop);
    }
  }


  fprintf(f,"\nBeginRange2\n");  
  printf("\nBeginRange2\n");  

  step=(smax2-smin2)/(NS2-1);
  for(i=0;i<NS2;i++)
  {
    sloop=(smin2+i*step);
    printf("%.5f \n",sloop);
    for(j=0;j<NCOS;j++)
    {
      costhetloop=-1.+1.0/NCOS+j*2./(NCOS);
 
      iew=0;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);
      borndivv=divv;

      iew=1;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);

      for(k=0;k<4;k++)
	for(l=0;l<4;l++)
	  fprintf(f,"%.5f ",R[k][l]);
      fprintf(f,"%.5f ",divv);
      fprintf(f,"%.5f \n",borndivv);
    }
  }


  fprintf(f,"\nBeginRange3\n");  
  printf("\nBeginRange3\n");  

  step=(smax3-smin3)/(NS3-1);
  for(i=0;i<NS3;i++)
  {
    sloop=(smin3+i*step);
    printf("%.5f \n",sloop);
    for(j=0;j<NCOS;j++)
    {
      costhetloop=-1.+1.0/NCOS+j*2./(NCOS);

      iew=0;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);
      borndivv=divv;

      iew=1;
      flagset_(&iqed,&iew,&iborn,&gfscheme,&ifgg,&nc,&fc,&tlmu2);
      Rcalc(flav,sloop,costhetloop);

      for(k=0;k<4;k++)
	for(l=0;l<4;l++)
	  fprintf(f,"%.5f ",R[k][l]);
      fprintf(f,"%.5f ",divv);
      fprintf(f,"%.5f \n",borndivv);
    }
  }
  fprintf(f,"End\n");
  fclose(f);
  return 0;
}