#include <stdio.h>
/* This programme calculates the action density, norm-squared of the Weyl zero
   mode and the Polyakov loop of the SU(2) caloron at omega=w, rho=p. For b=0
   it computes the action density, using n mirrors to -approximately- take the
   wrapping around the boundaries into account. Usually we take n=0, but
   carefully tune the location of the constituents wrt the boundary to minimize
   wrapping around the boundary. For b=-1 the zero-mode density (anti-periodic
   in time) is listed. For b=-2 the zero-mode (periodic in time) is listed,
   obtained by replacing (nu1,r1) and (nu2,r2). For b<-2 these two zero modes
   are added. For b>0 (forcing n=0) the Polykov loop with a discretization of
   b steps in the time direction is listed. The data is listed in a format
   suitable for comparision with data on a ns^3Xnt lattice at fixed y (and t).*/

void exit();
double w,p,x1,y1,z1,x2,y2,z2,t,Pi=3.14159265358979324;
double A[4][4],R[4][4];
int ns,nt,n,b;
     
main(argc,argv)
int argc;
char **argv;
{
        double mir();
        double pol();
        double sqrt();
        double x,y,z,t,s;
        int mx,my,mz,mt;
        void message();

        if(argc != 13)message();
        if(sscanf(argv[1],"%le",&w)!= 1)message();
        if(sscanf(argv[2],"%le",&x1)!= 1)message();
        if(sscanf(argv[3],"%le",&y1)!= 1)message();
        if(sscanf(argv[4],"%le",&z1)!= 1)message();
        if(sscanf(argv[5],"%le",&x2)!= 1)message();
        if(sscanf(argv[6],"%le",&y2)!= 1)message();
        if(sscanf(argv[7],"%le",&z2)!= 1)message();
        if(sscanf(argv[8],"%le",&t)!= 1)message();
        if(sscanf(argv[9],"%d",&ns)!= 1)message();
        if(sscanf(argv[10],"%d",&nt)!= 1)message();
        if(sscanf(argv[11],"%d",&n)!= 1)message();
        if(sscanf(argv[12],"%d",&b)!= 1)message();

        s=1.0/nt; 
        x=y=z=0.0;
        p=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)); 

        if(b>0){
          n=0;
          t=0.0;

       /* Below we define the rotation that brings the two constituents 
          parallel to the z-axis, with the new z coordinate of the first 
          consituent larger than that of the second. It is in this form 
          that the analytic formula for A used for computing the Polyakov
          loop are given!!! */

       /* vec u=(vec y2-vec y1)/|vec y1-vec y2|, vec R[2]=vec u wedge vec e1
          (R[2][1]=|vec R[2]|, using the 1st component of vec R[2] vanishes).
          and vec R[2] wedge vec u form the rows of the rotation matrix to put
          the constituents parallel to the z-axis. */

          R[3][1]=(x1-x2)/p;
          R[3][2]=(y1-y2)/p;
          R[3][3]=(z1-z2)/p;

          R[2][2]=z1-z2;
          R[2][3]=y2-y1;
          R[2][1]=sqrt(R[2][2]*R[2][2]+R[2][3]*R[2][3]);
          if(R[2][1]/p<1.0e-10){
              fprintf(stderr,"interchange x and z axis\n");
              exit();
          }
          R[2][2] /=R[2][1];
          R[2][3] /=R[2][1];
  
          /* we are now free to use R[2][1] for the norm of vec u wedge vec R[2]
             as long as we put it back to zero at the end */

          R[1][1]=R[2][2]*R[3][3]-R[2][3]*R[3][2];
          R[1][2]=R[2][3]*R[3][1];
          R[1][3]=-R[2][2]*R[3][1];
          R[2][1]=sqrt(R[1][1]*R[1][1]+R[1][2]*R[1][2]+R[1][3]*R[1][3]);
          R[1][1] /=R[2][1];
          R[1][2] /=R[2][1];
          R[1][3] /=R[2][1];
  
          R[2][1]  =0.0;
        }

        printf("(* omega=%lg, rho=%lg, number of mirrors=%d\n",w,sqrt(p/Pi),n);
        printf("   x1,y1,z1,ns= %lg, %lg, %lg, %d\n",x1,y1,z1,ns); 
        printf("   x2,y2,z2,nt= %lg, %lg, %lg, %d\n\n",x2,y2,z2,nt); 
        printf("   Slice through the constituents, with y=0 t0=%lg\n",t);
        printf("   from x=z=0 to ns/nt in steps of 1/nt\n");
        printf("   Pick xi and yi to optimize distance to boundaries\n");
        printf("   in absence of mirror contributions. *)\n\n");

        printf("Profile[%d]={\n",b); 

        if(b<1){
          for(mx=0;mx<ns-1;mx++){
            printf("{\n"); 
            for(mz=0;mz<ns-1;mz++){
              printf("%18.12E,\n",mir(x+s*mx,y,z+s*mz,t));
            }
            printf("%18.12E},\n",mir(x+s*mx,y,z+s*(ns-1),t));
          }
          printf("{\n"); 
          for(mz=0;mz<ns-1;mz++){
            printf("%18.12E,\n",mir(x+s*(ns-1),y,z+s*mz,t));
          }
          printf("%18.12E}\n",mir(x+s*(ns-1),y,z+s*(ns-1),t));
          printf("};\n"); 
        } else {
          for(mx=0;mx<ns-1;mx++){
            printf("{\n"); 
            for(mz=0;mz<ns-1;mz++){
              printf("%lg,\n",pol(x+s*mx,y,z+s*mz,b));
            }
            printf("%lg},\n",pol(x+s*mx,y,z+s*(ns-1),b));
         }
         printf("{\n"); 
         for(mz=0;mz<ns-1;mz++){
           printf("%lg,\n",pol(x+s*(ns-1),y,z+s*mz,b));
         }
         printf("%lg}\n",pol(x+s*(ns-1),y,z+s*(ns-1),b));
         printf("};\n"); 
      }

      return(0);
}
     
void message(){
      fprintf(stderr,"usage: SU2-lat omega x1 y1 z1 x2 y2 z2 -t0 ns nt n b\n");
      exit();
}

double pol(x,y,z,n)
double x,y,z;
int n;
{
       void AF();
       double sqrt(),cos(),sin();
       double A1,A2,A3,Ar,U0,U1,U2,U3,V0,V1,V2,V3,W0,W1,W2,W3;
       double s,t,q,xp,yp,zp;
       int m;

       q=ns*(1.0/nt);
       /* We first need to rotate the two constituents to be parallel to
          the z-axis, with the new z coordinate of the first consituent
          larger than that of the second. The Polyakov loop is simply
          compute with a finite-element approximation. */
       xp=R[1][1]*(x-x2)+R[1][2]*(y-y2)+R[1][3]*(z-z2);
       yp=R[2][1]*(x-x2)+R[2][2]*(y-y2)+R[2][3]*(z-z2);
       zp=R[3][1]*(x-x2)+R[3][2]*(y-y2)+R[3][3]*(z-z2);

       t=0.00001;
       s=1.0/n;
       V0=1.0;V1=V2=V3=0.0;
       for(m=0;m<n;m++){
          AF(2*w,0.,0.,0.,0.,p,0.,xp,yp,zp,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;
       }
       return(V0);
}

double mir(x,y,z,t)
double x,y,z,t;
{
       double dAdAlnpsi();
       double Nzmsq();
       double q,sum;
       int m1,m2,m3;

       q=ns*(1.0/nt);
       sum=0.0;
       if(b==0){
           for(m1=-n;m1<n+1;m1++) for(m2=-n;m2<n+1;m2++) for(m3=-n;m3<n+1;m3++)
              sum +=-dAdAlnpsi(2*w,x1,x2,y1,y2,z1,z2,x+m1*q,y+m2*q,z+m3*q,t)/2;
       }
       if(b==-1){
           for(m1=-n;m1<n+1;m1++) for(m2=-n;m2<n+1;m2++) for(m3=-n;m3<n+1;m3++)
              sum +=Nzmsq(2*w,x1,x2,y1,y2,z1,z2,x+m1*q,y+m2*q,z+m3*q,t);
       }
       if(b==-2){
           for(m1=-n;m1<n+1;m1++) for(m2=-n;m2<n+1;m2++) for(m3=-n;m3<n+1;m3++)
              sum +=Nzmsq(1-2*w,x2,x1,y2,y1,z2,z1,x+m1*q,y+m2*q,z+m3*q,t);
       }
       if(b<-2){
           for(m1=-n;m1<n+1;m1++) for(m2=-n;m2<n+1;m2++) for(m3=-n;m3<n+1;m3++){
              sum +=Nzmsq(2*w,x1,x2,y1,y2,z1,z2,x+m1*q,y+m2*q,z+m3*q,t);
              sum +=Nzmsq(1-2*w,x2,x1,y2,y1,z2,z1,x+m1*q,y+m2*q,z+m3*q,t);
           }
       }
       return(sum);
}