#include <stdlib.h>

#include <string.h>

#include <stdio.h>

#include <math.h>







struct vec4  {

	double x,y,z,t;

	vec4 operator - (vec4 v)   {vec4 uu;uu.x=x-v.x;uu.y=y-v.y;uu.z=z-v.z;uu.t=t-v.t;return uu;};

	vec4 operator / (double s) {vec4 uu;uu.x=x/s;uu.y=y/s;uu.z=z/s;uu.t=t/s;return uu;};



};

#define frand(u) (u*(double)(rand()%10000)/10000.0)







vec4 fncA(double t,double x,double y,double z)

{

	vec4 v4;

	double k=M_PI/180.0;

	double phase=z*k - t*k; 

	v4.x=sin(phase)*100.0;

//v4.x=0.0;	

	v4.y=0.0;

	v4.z= cos(phase)*100.0;//a hullam mozgas iranya, mivel phase=z*k

	v4.t=-cos(phase)*100.0;//time



	return v4;

}

void testtensor()

{

	double x,y,z,t;

	vec4 dAx,dAy,dAz,dAw;

	double dt=0.01, dx=0.01, dy=0.01, dz=0.01;

	x=frand(1000.0);

	y=frand(1000.0);

	z=frand(1000.0);

	t=frand(1000.0);



	

//four gradient

	dAw=(fncA(t+dt,x   ,y   ,z   ) - fncA(t-dt,x   ,y   ,z   ))/(dt*2.0);

	dAx=(fncA(t   ,x+dx,y   ,z   ) - fncA(t   ,x-dx,y   ,z   ))/(dx*2.0);

	dAy=(fncA(t   ,x   ,y+dy,z   ) - fncA(t   ,x   ,y-dy,z   ))/(dy*2.0);

	dAz=(fncA(t   ,x   ,y   ,z+dz) - fncA(t   ,x   ,y   ,z-dz))/(dz*2.0);





	double tensorF[4][4];

//exterior derivative  Aij - Aji

	tensorF[0][0]=dAw.t - dAw.t;

	tensorF[0][1]=dAw.x - dAx.t;//Ex   

	tensorF[0][2]=dAw.y - dAy.t;//Ey

	tensorF[0][3]=dAw.z - dAz.t;//Ez



	tensorF[1][0]=dAx.t - dAw.x;

	tensorF[1][1]=dAx.x - dAx.x;

	tensorF[1][2]=dAx.y - dAy.x;//Bz

	tensorF[1][3]=dAx.z - dAz.x;



	tensorF[2][0]=dAy.t - dAw.y;

	tensorF[2][1]=dAy.x - dAx.y;

	tensorF[2][2]=dAy.y - dAy.y;

	tensorF[2][3]=dAy.z - dAz.y;//Bx



	tensorF[3][0]=dAz.t - dAw.z;

	tensorF[3][1]=dAz.x - dAx.z;//By  

	tensorF[3][2]=dAz.y - dAy.z;

	tensorF[3][3]=dAz.z - dAz.z;



	if(fabs(tensorF[0][3])>1e-10 || fabs(tensorF[1][2])>1e-10)

	{

		printf("E %.3f %.3f %e \n",tensorF[0][1],tensorF[0][2],tensorF[0][3]);

		printf("B %.3f %.3f %e \n",tensorF[2][3],tensorF[3][1],tensorF[1][2]);

		printf("\n");

	}

	else printf("OK\n");

	

	return;

}







int main()

{

	for(int i=0;i<100;i++) testtensor();



	return 0;

}