#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 t 
   and z. Constituents are placed at x=y=0, z=(1-2*w)*p*p*Pi and -2*w*p*p*Pi.*/

main(argc,argv)
int argc;
char **argv;
{
        double dAdAlnpsi(),Nzmsq(); 
        double w,p,x,y,z,t,s,st,nu1,z1,z2;
        double Pi=3.14159265358979324;
        int mt,mz,nt,nz,nn=1;
        void message();

        if(argc != 10)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",&s)!= 1)message();
        if(sscanf(argv[7],"%d",&nz)!= 1)message();
        if(sscanf(argv[8],"%d",&nt)!= 1)message();
        if(sscanf(argv[9],"%d",&nn)!= 1)message(); 
        t=-0.5;
        st=nn*1.0/nt;
        nu1=2*w;
        z2=-nu1*p*p*Pi;
        z1=(1-nu1)*p*p*Pi;

        printf("(* t-z profile for: omega,rho,x,y= %lg, %lg, %lg, %lg\n",
                w,p,x,y); /* nn is the number of time periods */
        printf("   z[1],step,nz,nt,nn= %lg, %lg, %d, %d, %d *)\n",z,s,nz,nt,nn);
        printf("Profile[0]={\n"); 
        for(mt=0;mt<nt;mt++){
           printf("{\n"); 
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     -dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*mz,t+mt*st)/2); 
           printf("%16.10E},\n",
                  -dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*nz,t+st*mt)/2); 
           }
           printf("{\n"); 
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     -dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*mz,t+st*nt)/2); 
           printf("%16.10E}\n",
                  -dAdAlnpsi(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*nz,t+st*nt)/2); 
        printf("};\n\n"); 

        printf("Profile[-1]={\n");
        for(mt=0;mt<nt;mt++){
           printf("{\n");
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*mz,t+mt*st));
           printf("%16.10E},\n",
                  Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*nz,t+st*mt));
           }
           printf("{\n");
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*mz,t+st*nt));
           printf("%16.10E}\n",
                  Nzmsq(nu1,0.,0.,0.,0.,z1,z2,x,y,z+s*nz,t+st*nt));
        printf("};\n\n");

        printf("Profile[-2]={\n");
        for(mt=0;mt<nt;mt++){
           printf("{\n");
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x,y,z+s*mz,t+mt*st));
           printf("%16.10E},\n",
                  Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x,y,z+s*nz,t+st*mt));
           }
           printf("{\n");
           for(mz=0;mz<nz;mz++)
              printf("%16.10E,\n",
                     Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x,y,z+s*mz,t+st*nt));
           printf("%16.10E}\n",
                  Nzmsq(1-nu1,0.,0.,0.,0.,z2,z1,x,y,z+s*nz,t+st*nt));
        printf("};\n");

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