/**********************************************************
 *                                                        *
 *       Compute force volume in Vlasov equation for      *
 *   a few specific points in 6 coil volume.              *
 *                                                        *
 *              Author = Mike Rosing                      *
 *               date  = Oct. 17, 2007                    *
 *                                                        *
 *********************************************************/

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

main()
{
    FILE *magdat, *efld, *forcedat;
    double bx, by, bz, ux, uy, uz, fx, fy, fz;
    double ex, ey, ez, volts, amps, me, ce, param;
    int i, j, k, offset;
    char st[2];

    me = 9.10938189e-31;    // mass electron
    ce = 1.602176462e-19;   // charge of electron

    magdat = fopen("dipole.dat", "r");
    if(!magdat)
    {
	printf("can't open mag data\n");
	exit(0);
    }
    efld = fopen("e_field.dat", "r");
    if(!magdat)
    {
	printf("can't open e_field data\n");
	exit(0);
    }
    st[0] = 'n';
    while (st[0] == 'n')
    {
	printf("input i offset: ");
	scanf("%d", &i);
	printf("input j offset: ");
	scanf("%d", &j);
	printf("input k offset: ");
	scanf("%d", &k);
	offset = 3*sizeof(double)*(i + 256*(j + 256*k));
//	printf("desired set point is %d\n", offset);
	rewind(magdat);
	fseek( magdat, offset, SEEK_SET);
	rewind(efld);
	fseek( efld, offset, SEEK_SET);
//	offset = ftell(magdat);
//	printf("actuall set point is %d\n", offset);
	fread( (double*)&bx, sizeof(double), 1, magdat);
	fread( (double*)&by, sizeof(double), 1, magdat);
	fread( (double*)&bz, sizeof(double), 1, magdat);
	fread( (double*)&ex, sizeof(double), 1, efld);
	fread( (double*)&ey, sizeof(double), 1, efld);
	fread( (double*)&ez, sizeof(double), 1, efld);

	printf(" mag field values are bx = %lg  by = %lg  bz = %lg\n", bx, by, bz);
	printf(" e field values are ex = %lg  ey = %lg  ez = %lg\n", ex, ey, ez);
	printf("hit n to retry, y to continue: ");
	scanf("%s", st);
    }

/*  now compute u vector space and u X B  */

    forcedat = fopen("force.dat", "w");
    if (! forcedat)
    {
	printf("can't create force file\n");
	exit(0);
    }
    printf("Voltage on grid: ");
    scanf ("%lf", &volts);
    printf("Amp turns in grid: ");
    scanf("%lf", &amps);
    param = sqrt(2.0*me*volts/ce)/amps/(2.0e-7);
    printf("containment parameter = %lf\n", param);

    for(k=0; k<128; k++)
    {
	uz = k/32.0 - 2.0;
	for(j=0; j<128; j++)
	{
	    uy = j/32.0 - 2.0;
	    for(i=0; i<128; i++)
	    {
		ux = i/32.0 - 2.0;
		fx = param*ex + uy*bz - uz*by;
		fy = param*ey + uz*bx - ux*bz;
		fz = param*ez + ux*by - uy*bx;
		fwrite( &fx, sizeof(double), 1, forcedat);
		fwrite( &fy, sizeof(double), 1, forcedat);
		fwrite( &fz, sizeof(double), 1, forcedat);
	    }
	}
    }
    fclose(forcedat);
}