#include <stdio.h>
/* This programme calculates the action density profile (Profile[0]) and the
   zero-mode densities (anti-periodic in time Profile[-1]; periodic in time
   Profile[-2]) and of the SU(2) caloron at omega=w, rho=p as a function of
   r=sqrt(x^2+y^2) and z. Constituents are placed at r=0, z=(1-2*w)*p*p*Pi 
   and -2*w*p*p*Pi. */

main(argc,argv)
int argc;
char **argv;
{
        double dAdAlnpsi(),Nzmsq(),jac();
        double w,p,x,y,z,t,s,sum,tt,nu1,z1,z2;
        double Pi=3.14159265358979324;
        int mx,mz,nx,nz;
        void message();

        if(argc != 9)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",&z)!= 1)message();
        if(sscanf(argv[5],"%le",&t)!= 1)message();
        if(sscanf(argv[6],"%le",&s)!= 1)message();
        if(sscanf(argv[7],"%d",&nx)!= 1)message();
        if(sscanf(argv[8],"%d",&nz)!= 1)message();
        y=0.0;
        sum=0.0;
        nu1=2*w;
        z2=-nu1*p*p*Pi;
        z1=(1-nu1)*p*p*Pi;

        printf("(* r-z caloron profile for: omega= %lg and rho=%lg,\n",w,p); 
        printf("   r[1],z[1],t,step,nx,nz= %lg, %lg, %lg, %lg, %d, %d *)\n",
                x,z,t,s,nx,nz); 
        printf("Profile[0]={\n"); 
        for(mx=0;mx<nx-1;mx++){
           printf("{\n"); 
           for(mz=0;mz<nz-1;mz++){
               tt=-dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x+s*mx,y,z+s*mz,t)/2; 
               sum +=jac(x+s*mx)*tt;
               printf("%16.10E,\n",tt); 
           }
           tt=-dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x+s*mx,y,z+s*(nz-1),t)/2; 
           sum +=jac(x+s*mx)*tt;
           printf("%16.10E},\n",tt);
           }
           printf("{\n"); 
           for(mz=0;mz<nz-1;mz++){
              tt=-dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x+s*(nx-1),y,z+s*mz,t)/2; 
              sum +=jac(x+s*(nx-1))*tt;
              printf("%16.10E,\n",tt);
              }
           tt=-dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x+s*(nx-1),y,z+s*(nz-1),t)/2; 
           sum +=jac(x+s*(nx-1))*tt;
           printf("%16.10E}\n",tt);
        printf("};\n"); 
        printf("St[%lg]=%16.10E*(8Pi^2);\n\n",t,sum*s*s/8/Pi); 

        sum=0.0;
        printf("Profile[-1]={\n");
        for(mx=0;mx<nx-1;mx++){
           printf("{\n");
           for(mz=0;mz<nz-1;mz++){
               tt=Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x+s*mx,y,z+s*mz,t);
               sum +=jac(x+s*mx)*tt;
               printf("%16.10E,\n",tt);
           }
           tt=Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x+s*mx,y,z+s*(nz-1),t);
           sum +=jac(x+s*mx)*tt;
           printf("%16.10E},\n",tt);
           }
           printf("{\n");
           for(mz=0;mz<nz-1;mz++){
              tt=Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x+s*(nx-1),y,z+s*mz,t);
              sum +=jac(x+s*(nx-1))*tt;
              printf("%16.10E,\n",tt);
              }
           tt=Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x+s*(nx-1),y,z+s*(nz-1),t);
           sum +=jac(x+s*(nx-1))*tt;
           printf("%16.10E}\n",tt);
        printf("};\n");
        printf("Nzmsq[%lg]=%16.10E;\n\n",t,sum*s*s*Pi);

        sum=0.0;
        printf("Profile[-2]={\n");
        for(mx=0;mx<nx-1;mx++){
           printf("{\n");
           for(mz=0;mz<nz-1;mz++){
               tt=Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x+s*mx,y,z+s*mz,t);
               sum +=jac(x+s*mx)*tt;
               printf("%16.10E,\n",tt);
           }
           tt=Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x+s*mx,y,z+s*(nz-1),t);
           sum +=jac(x+s*mx)*tt;
           printf("%16.10E},\n",tt);
           }
           printf("{\n");
           for(mz=0;mz<nz-1;mz++){
              tt=Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x+s*(nx-1),y,z+s*mz,t);
              sum +=jac(x+s*(nx-1))*tt;
              printf("%16.10E,\n",tt);
              }
           tt=Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x+s*(nx-1),y,z+s*(nz-1),t);
           sum +=jac(x+s*(nx-1))*tt;
           printf("%16.10E}\n",tt);
        printf("};\n");
        printf("Nzmsq[%lg]=%16.10E;\n",t,sum*s*s*Pi);

        return(0);
}
     
void message(){
        fprintf(stderr,"usage: SU2-r-z omega rho r z t s nr nz\n");
        exit();
}

double jac(x)
double x;
{
         double sqrt();

         return(sqrt(x*x));
}