#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<fftlib.h>


#define VN 4       
#define pi M_PI  
#define M 4 
#define DN 512 
#define delta 0.2
#define n 11 


#define theta pi/4.0
#define shift pi/4.0
#define base pi/4.0
#define epsilon 0.0


void uv_p(double *g, double *vx, double *vy, double *vz, double *uvx, double *uvy, double *uvz ); 
void uv_s(double *g, double *sx[], double *sy[], double *sz[], double *vx, double *vy, double *vz, double *usx[], double *usy[], double *usz[] ); 

void us_ps(double *w, double *vx, double*vy, double *vz, double *sx[], double *sy[], double *sz[], double *uvx1, double *uvy1, double *uvz1, double *usx1[], double *usy1[], double *usz1[]);


void smallx(double *w, double *g, double *vx, double *vy, double *vz, double *fx, double *fy, double *fz, double *uvx, double *uvy, double *uvz, double *uvx1, double *uvy1, double *uvz1, double *sx[], double *sy[], double *sz[], double *usx1[], double *usy1[], double *usz1[]); 

void largex(double *w, double *g, double *sx[], double *sy[], double *sz[], double *gx[], double *gy[], double *gz[], double *usx[], double *usy[], double *usz[], double *usx1[], double *usy1[], double *usz1[], double *vx, double *vy, double *vz, double *uvx1, double *uvy1, double *uvz1); 


void rk4_m(double a, double *w, double *vx, double *vy, double *vz, double *uvx, double *uvy, double *uvz, double *uvx1, double *uvy1, double *uvz1, double *usx[], double *usy[], double *usz[], double *usx1[], double *usy1[], double *usz1[], double *sx[], double *sy[], double *sz[], double *g);

int main(){

int i,k;

double a=0.0;


double *vx, *vy, *vz, *g;
double z[M],*w;
double **sx, **sy, **sz; 
double *uvx, *uvy, *uvz;
double **usx, **usy, **usz;
double *uvx1, *uvy1, *uvz1;
double **usx1, **usy1, **usz1;
double **gx, **gy, **gz;
double *fx, *fy, *fz;


vx=FFTmallocDbl1D(VN);
vy=FFTmallocDbl1D(VN);
vz=FFTmallocDbl1D(VN);
g=FFTmallocDbl1D(VN);
w=FFTmallocDbl1D(M);
uvx=FFTmallocDbl1D(VN);
uvy=FFTmallocDbl1D(VN);
uvz=FFTmallocDbl1D(VN);
uvx1=FFTmallocDbl1D(VN);
uvy1=FFTmallocDbl1D(VN);
uvz1=FFTmallocDbl1D(VN);
fx=FFTmallocDbl1D(VN);
fy=FFTmallocDbl1D(VN);
fz=FFTmallocDbl1D(VN);



sx=FFTmallocDbl2D(M,DN);
sy=FFTmallocDbl2D(M,DN);
sz=FFTmallocDbl2D(M,DN);
usx=FFTmallocDbl2D(M,DN);
usy=FFTmallocDbl2D(M,DN);
usz=FFTmallocDbl2D(M,DN);
usx1=FFTmallocDbl2D(M,DN);
usy1=FFTmallocDbl2D(M,DN);
usz1=FFTmallocDbl2D(M,DN);
gx=FFTmallocDbl2D(M,DN);
gy=FFTmallocDbl2D(M,DN);
gz=FFTmallocDbl2D(M,DN);



for(i=0; i<M; i++){
   //z[i]=1.0-(2.0*(double)(i+1)/(M+1));
   w[i]=2.0/(double)(M+1);
  }


for(i=0;i<VN;i++){
        vz[i]=cos(base)+epsilon*sin(2.0*pi*(double)i/(double)VN);
	vx[i]=sqrt(1.0-vz[i]*vz[i]) * cos(2.0*pi*(double)i/(double)VN+ shift );
  	vy[i]=sqrt(1.0-vz[i]*vz[i]) * sin(2.0*pi*(double)i/(double)VN+ shift );
    	g[i]=4.0*sqrt(2.0)/3.0*pi;
	uvx[i]=0.0;
	uvy[i]=0.0;
	uvz[i]=0.0;
	uvx1[i]=0.0;
	uvy1[i]=0.0;
	uvz1[i]=0.0;
	fx[i]=0.0;
	fy[i]=0.0;
	fz[i]=0.0;
	printf("%lf %lf %lf\n",vx[i],vy[i],vz[i]);
}



for (i=0;i<M;i++){     
	for(k=0;k<DN;k++){
	         sz[i][k] = 1.0-2.0*(double)(i+1)/(double)(M+1);
                 sx[i][k] = sqrt(1.0-sz[i][k]*sz[i][k])*cos(2.0*pi*(double)k/(double)DN);
                 sy[i][k] = sqrt(1.0-sz[i][k]*sz[i][k])*sin(2.0*pi*(double)k/(double)DN);
	 	usx[i][k]=0.0;
                usy[i][k]=0.0;
		usz[i][k]=0.0;
		usx1[i][k]=0.0;
		usy1[i][k]=0.0;
		usz1[i][k]=0.0;
		gx[i][k]=0.0;
		gy[i][k]=0.0;
		gz[i][k]=0.0;
		printf("%lf %lf %lf\n",sx[i][k],sy[i][k],sz[i][k]);
		}
}


 
rk4_m(a, w, vx, vy, vz, uvx, uvy, uvz, uvx1, uvy1, uvz1, usx, usy, usz, usx1, usy1, usz1, sx, sy, sz, g);



FFTfreeDbl1D(vx);
FFTfreeDbl1D(vy);
FFTfreeDbl1D(vz);
FFTfreeDbl1D(g);
FFTfreeDbl1D(w);
FFTfreeDbl1D(uvx);
FFTfreeDbl1D(uvy);
FFTfreeDbl1D(uvz);
FFTfreeDbl1D(uvx1);
FFTfreeDbl1D(uvy1);
FFTfreeDbl1D(uvz1);
FFTfreeDbl1D(fx);
FFTfreeDbl1D(fy);
FFTfreeDbl1D(fz);



FFTfreeDbl2D(sx);
FFTfreeDbl2D(sy);
FFTfreeDbl2D(sz);
FFTfreeDbl2D(usx);
FFTfreeDbl2D(usy);
FFTfreeDbl2D(usz);
FFTfreeDbl2D(usx1);
FFTfreeDbl2D(usy1);
FFTfreeDbl2D(usz1);
FFTfreeDbl2D(gx);
FFTfreeDbl2D(gy);
FFTfreeDbl2D(gz);


return 0;

}


void uv_p(double *g, double *vx, double *vy, double *vz, double *uvx, double *uvy, double *uvz ){

int i,j;

for(i=0;i<VN;i++){
	uvx[i]=0.0;
	uvy[i]=0.0;
	uvz[i]=0.0;
	}

for(i=0; i<VN; i++){
     for(j=0; j<VN; j++){
		if(i!=j){
	    uvx[i] += g[j]*(vy[j]*vz[i]-vz[j]*vy[i])/(1.0-(vx[j]*vx[i]+vy[j]*vy[i]+vz[j]*vz[i]))/(4.0*pi);   
            uvy[i] += g[j]*(vz[j]*vx[i]-vx[j]*vz[i])/(1.0-(vx[j]*vx[i]+vy[j]*vy[i]+vz[j]*vz[i]))/(4.0*pi); 
            uvz[i] += g[j]*(vx[j]*vy[i]-vy[j]*vx[i])/(1.0-(vx[j]*vx[i]+vy[j]*vy[i]+vz[j]*vz[i]))/(4.0*pi); 
	    	}
	   }
}

}



void uv_s(double *g, double *sx[], double *sy[], double *sz[], double *vx, double *vy, double *vz, double *usx[], double *usy[], double *usz[] ){

int i,j,k;

for(i=0;i<M; i++){
	for(k=0;k<DN;k++){
		usx[i][k]=0.0;
		usy[i][k]=0.0;
		usz[i][k]=0.0;
	}
}


for(i=0;i<M;i++){
  	for(k=0;k<DN;k++){
	  for(j=0;j<VN;j++){	  
        usx[i][k] += g[j]*(vy[j]*sz[i][k]-vz[j]*sy[i][k])/(1.0+(delta*delta)-(vx[j]*sx[i][k]+vy[j]*sy[i][k]+vz[j]*sz[i][k]))/(4.0*pi);
	usy[i][k] += g[j]*(vz[j]*sx[i][k]-vx[j]*sz[i][k])/(1.0+(delta*delta)-(vx[j]*sx[i][k]+vy[j]*sy[i][k]+vz[j]*sz[i][k]))/(4.0*pi);
	usz[i][k] += g[j]*(vx[j]*sy[i][k]-vy[j]*sx[i][k])/(1.0+(delta*delta)-(vx[j]*sx[i][k]+vy[j]*sy[i][k]+vz[j]*sz[i][k]))/(4.0*pi);
	  }
	}

}


}


void us_ps(double *w, double *vx, double*vy, double *vz, double *sx[], double *sy[], double *sz[], double *uvx1, double *uvy1, double *uvz1, double *usx1[], double *usy1[], double *usz1[]){ 

int i,j,k,l;
Complex **c1, **c2, **d1, **d2, **e1, **e2;
double R;

/*??????*/
c1 = FFTmallocComplex2D(M,DN); /* ?????????W*/
c2 = FFTmallocComplex2D(M,DN);

d1 = FFTmallocComplex2D(M,DN);/*?????????W*/
d2 = FFTmallocComplex2D(M,DN);

e1 = FFTmallocComplex2D(M,DN);/*????????W*/
e2 = FFTmallocComplex2D(M,DN);



//*?????????W*/
  for (i = 0; i < M; i++) {
    for(k = 0; k < DN; k++){
    c1[i][k].r = sx[i][k];
    c1[i][k].i = 0.0;
     }
  }

  
  /* FFT (????? ??FFT ??  ??w???????????\??
     ?e?L?X?g???`??-1 ?????? */
  for(i=0;i<M;i++){
	  FFT1D(DN, -1, c1[i]);
  	
	for(k=0; k<DN; k++){
		if(0 <= k && k < (DN/2)){
			 c2[i][k].r = 2.0*pi*k*(-1.0)*c1[i][k].i;
			 c2[i][k].i = 2.0*pi*k*c1[i][k].r;
                }else if((DN/2) <= k && k < DN){
                        c2[i][k].r = 2.0*pi*(k-DN)*(-1.0)*c1[i][k].i;
                        c2[i][k].i = 2.0*pi*(k-DN)*c1[i][k].r;
	         }
	 
       }
  
	 
  /* FFT ?t???*/
  	FFT1D(DN, 1, c2[i]);
  	for (k = 0; k < DN; k++) {
    		c2[i][k].r /= (DN);
  	}
   }



/*?????????W?????*/
for (i = 0; i < M; i++) {
    for(k = 0; k < DN; k++){
    d1[i][k].r = sy[i][k];
    d1[i][k].i = 0.0;
    }
  }

  /* FFT (????? ??FFT ??  ??w???????????\??
     ?e?L?X?g???`??-1 ?????? */
  for(i=0;i<M;i++){
	  FFT1D(DN, -1, d1[i]);
  	
	for(k=0; k<DN; k++){
		if(0 <= k && k < (DN/2)){
			 d2[i][k].r=2.0*pi*k*(-1.0)*d1[i][k].i;
			 d2[i][k].i=2.0*pi*k*d1[i][k].r;
                }else if((DN/2) <= k && k < DN){
                        d2[i][k].r=2.0*pi*(k-DN)*(-1.0)*d1[i][k].i;
                        d2[i][k].i=2.0*pi*(k-DN)*d1[i][k].r;
		 }
       }
  
 
  /* FFT ?t???*/
  	FFT1D(DN, 1, d2[i]);
  	for (k = 0; k < DN; k++) {
    		d2[i][k].r /= (DN);/*?????????W*/
  	}

 }


/*????????W*/
  for (i = 0; i < M; i++) {
    for(k = 0; k < DN; k++){
    e1[i][k].r = sz[i][k];
    e1[i][k].i = 0.0;
    }
  }

  /* FFT (????? ??FFT ??  ??w???????????\??
     ?e?L?X?g???`??-1 ?????? */
  for(i=0;i<M;i++){
	  FFT1D(DN, -1, e1[i]);
  	

	for(k=0; k<DN; k++){
		if(0 <= k && k < (DN/2)){
			 e2[i][k].r = 2.0*pi*k*(-1.0)*e1[i][k].i;
			 e2[i][k].i = 2.0*pi*k*e1[i][k].r;
                }else if((DN/2) <= k && k < DN){
                        e2[i][k].r=2.0*pi*(k-DN)*(-1.0)*e1[i][k].i;
                        e2[i][k].i=2.0*pi*(k-DN)*e1[i][k].r;
		 }
       }
  
 
  /* FFT ?t???*/
  	FFT1D(DN, 1, e2[i]);
  	for (k = 0; k < DN; k++) {
    		e2[i][k].r /= (DN);/*?????????W???????*/
  	}
   }


for(i=0;i<VN;i++){
	uvx1[i]=0.0;
	uvy1[i]=0.0;
	uvz1[i]=0.0;
}

for(i=0;i<M;i++){
	for(j=0;j<DN;j++){
		usx1[i][j]=0.0;
		usy1[i][j]=0.0;
		usz1[i][j]=0.0;
	}
	}


for(j=0;j<VN;j++){
   for(i=0;i<M;i++){
	for(k=1;k<DN;k++){
	  R=(vx[j]-sx[i][k])*(vx[j]-sx[i][k])+(vy[j]-sy[i][k])*(vy[j]-sy[i][k])+(vz[j]-sz[i][k])*(vz[j]-sz[i][k])+(delta*delta);
		uvx1[j] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * c2[i][k].r) ;
		uvy1[j] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * d2[i][k].r) ;
		uvz1[j] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * e2[i][k].r) ;
		}
	}
}


for(j=0;j<M;j++){
	for(l=0;l<DN;l++){	
		for(i=0;i<M;i++){  
  			for(k=1;k<DN;k++){
			  if(i==j && k==l)continue;
	R=(sx[j][l]-sx[i][k])*(sx[j][l]-sx[i][k])+(sy[j][l]-sy[i][k])*(sy[j][l]-sy[i][k])+(sz[j][l]-sz[i][k])*(sz[j][l]-sz[i][k]);
				usx1[j][l] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * c2[i][k].r ) ;
				usy1[j][l] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * d2[i][k].r ) ;
				usz1[j][l] += (-1.0)/(4.0*pi*DN) * w[i] * (log(R) * e2[i][k].r ) ;
					
				}
			}
	}
}




FFTfreeComplex2D(c1);
FFTfreeComplex2D(c2);

FFTfreeComplex2D(d1);
FFTfreeComplex2D(d2);

FFTfreeComplex2D(e1);
FFTfreeComplex2D(e2);




}


void smallx(double *w, double *g, double *vx, double *vy, double *vz, double *fx, double *fy, double *fz, double *uvx, double *uvy, double *uvz, double *uvx1, double *uvy1, double *uvz1, double *sx[], double *sy[], double *sz[], double *usx1[], double *usy1[], double *usz1[]){ 

int i;

for(i=0;i<VN;i++){
	fx[i] = 0.0;
	fy[i] = 0.0;
	fz[i] = 0.0;
}

uv_p(g, vx, vy, vz, uvx, uvy, uvz);
us_ps(w, vx, vy, vz, sx, sy, sz, uvx1, uvy1, uvz1, usx1, usy1, usz1);

/*for(i=0;i<VN;i++){
	printf("%lf %lf %lf %lf %lf %lf\n",uvx[i], uvy[i], uvz[i],uvx1[i], uvy1[i], uvz1[i]);
	}*/


for(i=0;i<VN;i++){
	fx[i] = uvx[i] + uvx1[i];
	fy[i] = uvy[i] + uvy1[i];
	fz[i] = uvz[i] + uvz1[i];

	/*printf("%lf %lf %lf\n",fx[i],fy[i],fz[i]);*/

}


}


void largex(double *w, double *g, double *sx[], double *sy[], double *sz[], double *gx[], double *gy[], double *gz[], double *usx[], double *usy[], double *usz[], double *usx1[], double *usy1[], double *usz1[], double *vx, double *vy, double *vz, double *uvx1, double *uvy1, double *uvz1){ 

int i,j;

for(i=0;i<M;i++){
	for(j=0;j<DN;j++){
		gx[i][j] = 0.0;
		gy[i][j] = 0.0;
		gz[i][j] = 0.0;
	}
}

uv_s(g, sx, sy, sz, vx, vy, vz, usx, usy, usz);
us_ps(w, vx, vy, vz, sx, sy, sz, uvx1, uvy1, uvz1, usx1, usy1, usz1);

/*for(i=0;i<M;i++){
  for(j=0;j<DN;j++){
	printf("%lf %lf %lf %lf %lf %lf\n",usx[i][j], usy[i][j], usz[i][j],usx1[i][j], usy1[i][j], usz1[i][j]);
}
}*/


for(i=0;i<M;i++){
	for(j=0;j<DN;j++){
		gx[i][j] = usx[i][j] + usx1[i][j];
		gy[i][j] = usy[i][j] + usy1[i][j];
		gz[i][j] = usz[i][j] + usz1[i][j];

		//	printf("%lf %lf %lf\n",gx[i][j],gy[i][j],gz[i][j]);
	}
	/*printf("\n");*/
}

}


void rk4_m(double a, double *w, double *vx, double *vy, double *vz, double *uvx, double *uvy, double *uvz, double *uvx1, double *uvy1, double *uvz1, double *usx[], double *usy[], double *usz[], double *usx1[], double *usy1[], double *usz1[], double *sx[], double *sy[], double *sz[], double *g){

double *k1x, *k2x, *k3x, *k4x, *tmpx;
double *k1y, *k2y, *k3y, *k4y, *tmpy;
double *k1z, *k2z, *k3z, *k4z, *tmpz;

double **l1x, **l1y, **l1z, **tmpsx;
double **l2x, **l2y, **l2z, **tmpsy;
double **l3x, **l3y, **l3z, **tmpsz;
double **l4x, **l4y, **l4z;

double t,h;

//int l=0;
int  i,j,k;



k1x = FFTmallocDbl1D(VN);   
k2x = FFTmallocDbl1D(VN);
k3x = FFTmallocDbl1D(VN);   
k4x = FFTmallocDbl1D(VN); 
tmpx = FFTmallocDbl1D(VN);

k1y = FFTmallocDbl1D(VN);   
k2y = FFTmallocDbl1D(VN); 
k3y = FFTmallocDbl1D(VN);   
k4y = FFTmallocDbl1D(VN);
tmpy = FFTmallocDbl1D(VN);

k1z = FFTmallocDbl1D(VN);   
k2z = FFTmallocDbl1D(VN); 
k3z = FFTmallocDbl1D(VN);   
k4z = FFTmallocDbl1D(VN); 
tmpz = FFTmallocDbl1D(VN);

l1x = FFTmallocDbl2D(M,DN);   
l2x = FFTmallocDbl2D(M,DN);
l3x = FFTmallocDbl2D(M,DN);   
l4x = FFTmallocDbl2D(M,DN); 
tmpsx = FFTmallocDbl2D(M,DN);

l1y = FFTmallocDbl2D(M,DN);   
l2y = FFTmallocDbl2D(M,DN); 
l3y = FFTmallocDbl2D(M,DN);   
l4y = FFTmallocDbl2D(M,DN);
tmpsy = FFTmallocDbl2D(M,DN);

l1z = FFTmallocDbl2D(M,DN);   
l2z = FFTmallocDbl2D(M,DN); 
l3z = FFTmallocDbl2D(M,DN);   
l4z = FFTmallocDbl2D(M,DN); 
tmpsz = FFTmallocDbl2D(M,DN);
  

 h = 0.05;            t = a;




for ( i = 0 ; i < n ; i++) {  

  /* k1 */
      	smallx(w, g, vx, vy, vz, k1x, k1y, k1z, uvx, uvy, uvz, uvx1, uvy1, uvz1, sx, sy, sz, usx1, usy1, usz1); 
	largex(w, g, sx, sy, sz, l1x, l1y, l1z, usx, usy, usz, usx1, usy1, usz1, vx, vy, vz, uvx1, uvy1, uvz1);

	

	/*for(j=0;j<VN;j++){
	  printf("%lf %lf %lf\n",k1x[j],k1y[j],k1z[j]);
	  }*/

	for( j = 0; j < VN; j++){
		 tmpx[j] = vx[j] + h*k1x[j]/2.0;
	 	 tmpy[j] = vy[j] + h*k1y[j]/2.0;
		 tmpz[j] = vz[j] + h*k1z[j]/2.0;
	}
	for( j = 0; j < M; j++){ 
		for( k = 0; k < DN; k++){
			tmpsx[j][k] = sx[j][k] + h*l1x[j][k]/2.0;
			tmpsy[j][k] = sy[j][k] + h*l1y[j][k]/2.0;
			tmpsz[j][k] = sz[j][k] + h*l1z[j][k]/2.0;
		}
	}

      /* k2 */
	smallx(w, g, tmpx, tmpy, tmpz, k2x, k2y, k2z, uvx, uvy, uvz, uvx1, uvy1, uvz1, tmpsx, tmpsy, tmpsz, usx1, usy1, usz1);
        largex(w, g, tmpsx, tmpsy, tmpsz, l2x, l2y, l2z, usx, usy, usz, usx1, usy1, usz1, tmpx, tmpy, tmpz, uvx1, uvy1, uvz1);

		
	
	/*for(j=0;j<VN;j++){
	  printf("%lf %lf %lf\n",k2x[j],k2y[j],k2z[j]);
	  }*/


	for( j = 0; j < VN; j++){
		 tmpx[j] = vx[j] + h*k2x[j]/2.0;
		 tmpy[j] = vy[j] + h*k2y[j]/2.0;
		 tmpz[j] = vz[j] + h*k2z[j]/2.0;
	}
	for( j = 0; j < M; j++){ 
		for( k = 0; k < DN; k++){
			tmpsx[j][k] = sx[j][k] + h*l2x[j][k]/2.0;
			tmpsy[j][k] = sy[j][k] + h*l2y[j][k]/2.0;
			tmpsz[j][k] = sz[j][k] + h*l2z[j][k]/2.0;
		}
	}

      /* k3*/
        smallx(w, g, tmpx, tmpy, tmpz, k3x, k3y, k3z, uvx, uvy, uvz, uvx1, uvy1, uvz1, tmpsx, tmpsy, tmpsz, usx1, usy1, usz1);
	largex(w, g, tmpsx, tmpsy, tmpsz, l3x, l3y, l3z, usx, usy, usz, usx1, usy1, usz1, tmpx, tmpy, tmpz, uvx1, uvy1, uvz1);
        

	/*for(j=0;j<VN;j++){
	  printf("%lf %lf %lf\n",k3x[j],k3y[j],k3z[j]);
	  }*/


	for( j = 0; j < VN; j++ ){
		 tmpx[j] = vx[j] + h*k3x[j];
	         tmpy[j] = vy[j] + h*k3y[j];
                 tmpz[j] = vz[j] + h*k3z[j];
	}
	for( j = 0; j < M; j++){ 
		for( k = 0; k < DN; k++){
			tmpsx[j][k] = sx[j][k] + h*l3x[j][k];
			tmpsy[j][k] = sy[j][k] + h*l3y[j][k];
			tmpsz[j][k] = sz[j][k] + h*l3z[j][k];
		}
	}

      /* k4*/
        smallx(w, g, tmpx, tmpy, tmpz, k4x, k4y, k4z, uvx, uvy, uvz, uvx1, uvy1, uvz1, tmpsx, tmpsy, tmpsz, usx1, usy1, usz1);
	largex(w, g, tmpsx, tmpsy, tmpsz, l4x, l4y, l4z, usx, usy, usz, usx1, usy1, usz1, tmpx, tmpy, tmpz, uvx1, uvy1, uvz1);
    
	

	/*for(j=0;j<VN;j++){
	  printf("%lf %lf %lf\n",k4x[j],k4y[j],k4z[j]);
	  }*/

	printf("\t%lf\n",t);	
	
	for( j = 0; j < VN; j++ ){
		vx[j] = vx[j] + h*( k1x[j] + 2.0*k2x[j] + 2.0*k3x[j] + k4x[j] )/6.0;
        	vy[j] = vy[j] + h*( k1y[j] + 2.0*k2y[j] + 2.0*k3y[j] + k4y[j] )/6.0;
       		vz[j] = vz[j] + h*( k1z[j] + 2.0*k2z[j] + 2.0*k3z[j] + k4z[j] )/6.0;
		printf("%lf %lf %lf\n",vx[j],vy[j],vz[j]);	 
	       }

	for( j = 0; j < M; j++){ 
		for( k = 0; k < DN; k++){
			sx[j][k] = sx[j][k] + h*(l1x[j][k] + 2.0*l2x[j][k] + 2.0*l3x[j][k] + l4x[j][k] )/6.0;
			sy[j][k] = sy[j][k] + h*(l1y[j][k] + 2.0*l2y[j][k] + 2.0*l3y[j][k] + l4y[j][k] )/6.0;
			sz[j][k] = sz[j][k] + h*(l1z[j][k] + 2.0*l2z[j][k] + 2.0*l3z[j][k] + l4z[j][k] )/6.0;
			printf("%lf %lf %lf\n",sx[j][k],sy[j][k],sz[j][k]);			}
	}


	/*if(i==l){
	   printf("%lf\n",t);
	   for(j=0;j<M;j++){
			for(k=0;k<DN;k++){
				printf("%lf %lf %lf\n",sx[j][k],sy[j][k],sz[j][k]);
			}
			printf("\n");
	   }
	   for(j=0;j<VN;j++){
			printf("%lf %lf %lf\n",vx[j],vy[j],vz[j]);
	    }
	l+=1;
	}*/

	t+=h;


}






FFTfreeDbl1D(k1x);
FFTfreeDbl1D(k2x);
FFTfreeDbl1D(k3x);
FFTfreeDbl1D(k4x);
FFTfreeDbl1D(tmpx);

FFTfreeDbl1D(k1y);
FFTfreeDbl1D(k2y);
FFTfreeDbl1D(k3y);
FFTfreeDbl1D(k4y);
FFTfreeDbl1D(tmpy);

FFTfreeDbl1D(k1z);
FFTfreeDbl1D(k2z);
FFTfreeDbl1D(k3z);
FFTfreeDbl1D(k4z);
FFTfreeDbl1D(tmpz);

FFTfreeDbl2D(l1x);
FFTfreeDbl2D(l2x);
FFTfreeDbl2D(l3x);
FFTfreeDbl2D(l4x);
FFTfreeDbl2D(tmpsx);

FFTfreeDbl2D(l1y);
FFTfreeDbl2D(l2y);
FFTfreeDbl2D(l3y);
FFTfreeDbl2D(l4y);
FFTfreeDbl2D(tmpsy);

FFTfreeDbl2D(l1z);
FFTfreeDbl2D(l2z);
FFTfreeDbl2D(l3z);
FFTfreeDbl2D(l4z);
FFTfreeDbl2D(tmpsz);



}