#include <stdio.h>
/* This programme calculates the field and action density of the SU(2) caloron
   with the two constituents on the line x=y=0, placed at z=yz1>0 and -yz1.
   It also calculates the zero-mode density (anti-periodic in time) and the 
   trace of the Polyakov loop. The loop is straight starting at t0 (=0.00001, 
   to avoid hitting the gauge singularity) in nt steps at (x,y,z). */

double A[4][4],F[4][4][4];

main(argc,argv)
int argc;
char **argv;
{
        void message();
        void AF();
        double Nzmsq();
        double dAdAlnpsi();
        double sqrt(),cos(),sin();
        double Pi=3.14159265358979324;
        double yx1,yx2,yy1,yy2,yz1,yz2,nu1,x,y,z,t,s,w,p;
        double A1,A2,A3,Ar,U0,U1,U2,U3,V0,V1,V2,V3,W0,W1,W2,W3;
        double B1sq=0.0,B2sq=0.0,B3sq=0.0;
        int a,n,m,nt;

        if(argc != 8)message();
        if(sscanf(argv[1],"%le",&w)!= 1)message();
        if(sscanf(argv[2],"%le",&p)!= 1)message();
        if(sscanf(argv[3],"%le",&x)!= 1)message();
        if(sscanf(argv[4],"%le",&y)!= 1)message();
        if(sscanf(argv[5],"%le",&z)!= 1)message();
        if(sscanf(argv[6],"%le",&t)!= 1)message();
        if(sscanf(argv[7],"%d",&nt)!= 1)message();
        if(w>0.5 || w<0.0){
          printf("omega=%lg is not between 0 and 1/2\n",w);
          exit();
        }
        yx1=yx2=yy1=yy2=0.0;
        nu1=2*w;
        yz1=(1-nu1)*p*p*Pi;
        yz2=-nu1*p*p*Pi;

     	printf("nu1=%lg,r1=(%lg,%lg,%lg),nu2=%lg,r2=(%lg,%lg,%lg),X=(x,y,z,t)=(%lg,%lg,%lg,%lg),nt=%d\n\n",
                nu1,yx1,yy1,yz1,1-nu1,yx2,yy2,yz2,x,y,z,t,nt);
        AF(nu1,yx1,yx2,yy1,yy2,yz1,yz2,x,y,z,t);
        for(a=1;a<4;a++) B1sq +=F[a][2][3]*F[a][2][3];
        for(a=1;a<4;a++) B2sq +=F[a][1][3]*F[a][1][3];
        for(a=1;a<4;a++) B3sq +=F[a][1][2]*F[a][1][2];
        for(n=0;n<4;n++){
     	printf("A[1][%d]=%13.6E, A[2][%d]=%13.6E, A[3][%d]=%13.6E,\n",
                n,A[1][n],n,A[2][n],n,A[3][n]);
        }
        printf("\n");
        for(a=1;a<4;a++){
     	printf("B[%d][1]=%13.6E, B[%d][2]=%13.6E, B[%d][3]=%13.6E\n",
                a,F[a][2][3],a,F[a][3][1],a,F[a][1][2]);
        }
        printf("\n");
     	printf("B1^2(X)=%12.6E, B2^2(X)=%12.6E, B3^2(X)=%12.6E\n\n",
                B1sq,B2sq,B3sq);
     	printf("SFsq(X)=%16.10E\n",B1sq+B2sq+B3sq);
     	printf("SPsi(X)=%16.10E\n\n",
               -dAdAlnpsi(nu1,yx1,yx2,yy1,yy2,yz1,yz2,x,y,z,t)/2);
        printf("Nzmsq(X)=%16.10E\n\n",
               Nzmsq(nu1,yx1,yx2,yy1,yy2,yz1,yz2,x,y,z,t));

        t=0.00001;
        s=1.0/nt;
        V0=1.0;V1=V2=V3=0.0;

        for(m=0;m<nt;m++){
            AF(nu1,yx1,yx2,yy1,yy2,yz1,yz2,x,y,z,t+m*s);
            A1=A[1][0];A2=A[2][0];A3=A[3][0];
            Ar=sqrt(A1*A1+A2*A2+A3*A3);
            U1=U2=U3=0.0;
            U0=sin(Ar*s/2);
            if(Ar!=0.){
                  U1=A1*U0/Ar;
                  U2=A2*U0/Ar;
                  U3=A3*U0/Ar;
            }
            U0=cos(Ar*s/2);
            W0=V0*U0-V1*U1-V2*U2-V3*U3;
            W1=V0*U1+V1*U0-V2*U3+V3*U2;
            W2=V0*U2+V2*U0-V3*U1+V1*U3;
            W3=V0*U3+V3*U0-V1*U2+V2*U1;
            V0=W0;V1=W1;V2=W2;V3=W3;
        }
     	printf("Pol(X)=%16.10E\n",V0);
        return(0);
}
     
void message(){
        fprintf(stderr,"usage: SU2 omega rho x y z t nt\n");
        exit();
}