/*************************************************************************
 *                                                                       *
 *    This program creates magnetic field data for plotting using line   *
 *  segments for currents.  The analytical formula is derived using      *
 *  a parametric form for the line segment and integrating the           *
 *  (J x r)/r^2 for each point along the line segment.  In principle the *
 *  entire array should be symmetric, but for simplicity I compute the   *
 *  full field.                                                          *
 *                                                                       *
 *                      Author = Mike Rosing                             *
 *                       date  = June 1, 2008                            *
 *                                                                       *
 ************************************************************************/


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

#define MAXDIM 256
//#define MAXDIM 16

typedef struct
{
    double x, y, z;
} VECTOR;

VECTOR Bfield[MAXDIM][MAXDIM][MAXDIM];

main()
{
    struct
    {
	int start, end;
    } list[12];
    
    VECTOR location, direction[12], endpoint[6];

    int i, ix, iy, iz, start;
    double x, y, z, xp, yp, zp, rp2, k, sq, t1, t2;
    FILE *saveb;

    endpoint[0].x = 0.0;
    endpoint[0].y = 0.0;
    endpoint[0].z = 1.0;
    endpoint[1].x = 1.0;
    endpoint[1].y = 0.0;
    endpoint[1].z = 0.0;
    endpoint[2].x = 0.0;
    endpoint[2].y = 1.0;
    endpoint[2].z = 0.0;
    endpoint[3].x = -1.0;
    endpoint[3].y = 0.0;
    endpoint[3].z = 0.0;
    endpoint[4].x = 0.0;
    endpoint[4].y = -1.0;
    endpoint[4].z = 0.0;
    endpoint[5].x = 0.0;
    endpoint[5].y = 0.0;
    endpoint[5].z = -1.0;
    list[0].start = 0;
    list[0].end = 1;
    list[1].start = 1;
    list[1].end = 4;
    list[2].start = 4;
    list[2].end = 5;
    list[3].start = 5;
    list[3].end = 1;
    list[4].start = 1;
    list[4].end = 2;
    list[5].start = 2;
    list[5].end = 0;
    list[6].start = 0;
    list[6].end = 3;
    list[7].start = 3;
    list[7].end = 2;
    list[8].start = 2;
    list[8].end = 5;
    list[9].start = 5;
    list[9].end = 3;
    list[10].start = 3;
    list[10].end = 4;
    list[11].start = 4;
    list[11].end = 0;

/*  generate direction vectors of line segments  */

    for(i=0; i<12; i++)
    {
	direction[i].x = endpoint[list[i].end].x - endpoint[list[i].start].x;
	direction[i].y = endpoint[list[i].end].y - endpoint[list[i].start].y;
	direction[i].z = endpoint[list[i].end].z - endpoint[list[i].start].z;
	printf("direction[%d] x = %f  y = %f  z = %f\n", i,direction[i].x,
	       direction[i].y, direction[i].z);
    }

/*  assume all line segments are same length.  Sweep over whole volume and 
    add field from each line segement to the volume element.
*/

    for(ix=0; ix<MAXDIM; ix++)
    {
	x = (ix - 256.0)/64.0;
//	x = (ix - 8.0)/4.0;
	    printf("x= %d\n", ix);
	for(iy=0; iy<MAXDIM; iy++)
	{
	    y = (iy - 128.0)/64.0;
//	    y = (iy - 8.0)/4.0;
//	    printf("  y = %d\n", iy);
	    for(iz=0; iz<MAXDIM; iz++)
	    {
		z = (iz - 128.0)/64.0;
//		z = (iz - 8.0)/4.0;
//		printf("       z = %d\n", iz);
		Bfield[ix][iy][iz].x = 0.0;
		Bfield[ix][iy][iz].y = 0.0;
		Bfield[ix][iy][iz].z = 0.0;
		for(i=0; i<12; i++)
//		i = 5;
		{
		    start = list[i].start;
		    xp = x - endpoint[start].x;
		    yp = y - endpoint[start].y;
		    zp = z - endpoint[start].z;
		    rp2 = xp*xp + yp*yp + zp*zp;
		    k = 2.0*(direction[i].x*xp + direction[i].y*yp + direction[i].z*zp);
//		    printf("k = %f  rp2 = %f \n", k, rp2);
		    sq = 8.0*rp2 - k*k;
		    if(sq < 1e-10) continue;
		    if (sq < 0.0)
		    {
			printf("at i = %d j=%d k=%d sq = %f\n", ix, iy, iz, sq);
			exit(0);
		    }
		    sq = sqrt(sq);
		    t1 = atan2((k + 4.0),sq) - atan2(k,sq);
		    if (isnan(t1))
		      printf("i=%d ix=%d  iy=%d  iz=%d\n", i, ix, iy, iz);
//		    printf("t1 = %f\n", t1);
		    t2 = t1/sqrt(2.0)/sq;
//		    if(iy==127 && iz==127) printf("t2=%f\n",t2);
		    Bfield[ix][iy][iz].x += (direction[i].y*zp - direction[i].z*yp)*t2;
		    Bfield[ix][iy][iz].y += (direction[i].z*xp - direction[i].x*zp)*t2;
		    Bfield[ix][iy][iz].z += (direction[i].x*yp - direction[i].y*xp)*t2;
/*		    if(ix==iz && iy == iz) printf("bx = %f  by = %f  bz = %f\n",
					       Bfield[ix][iy][iz].x, Bfield[ix][iy][iz].y,
					       Bfield[ix][iy][iz].z);
*/
		}
	    }
	}
    }

    saveb = fopen("tombo_field.dat", "w");
    if(!saveb)
    {
	printf("can't create tombo_field.dat\n");
	exit(0);
    }
/*    for(ix=0; ix<MAXDIM; ix++)
    {
	for(iy=0; iy<MAXDIM; iy++)
	{
	    for(iz=0; iz<MAXDIM; iz++) fprintf(saveb, "%d %f %f %f\n", 
					       (ix*MAXDIM + iy)*MAXDIM+ iz, 
					       Bfield[ix][iy][iz].x,
					    Bfield[ix][iy][iz].y, Bfield[ix][iy][iz].z);
	}
	}*/
    fwrite(Bfield, sizeof(VECTOR), MAXDIM*MAXDIM*MAXDIM, saveb);
    fclose(saveb);
    printf("all done\n");
}