/******************************************************************************
 *                                                                            *
 *    Read in magnetic potential file and create text file that gnuplot       *
 *  can read.                                                                 *
 *                                                                            *
 *                       Author = Mike Rosing                                 *
 *                        Date  = Feb 16, 2008                                *
 *                                                                            *
 *****************************************************************************/

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

#define MAXSTEPS 25

/*  MAXWALL is the ratio of outer containment sphere radius to L (center
    of sphere to center of coil distance)  */

#define MAXWALL  (3.0)

/*  COILRADIUS is the ratio of coil radius to L  */

#define COILRADIUS (0.9)
#define NUMELMNTS (3*(MAXSTEPS)*(MAXSTEPS+1)*(MAXSTEPS+2)/6)

/*  Create global buffers for data storage so integrals are easy */

double MagPot[NUMELMNTS];

main()
{
    FILE  *readmag, *writedat;
    int numread;
    char filename[64], input[8];
    int i, j, k, idex, kdex;
    double across, r, phi, theta;
    double pi4;

    pi4 = M_PI/4.0;          // ubiquitous constant

/*  load grid potential from previous calculation  */

    printf("File to convert for gnuplot: ");
    scanf("%s", filename);
//    sprintf(filename, "electron_potential.dat");
    readmag = fopen(filename, "rb");
    if( !readmag)
    {
	printf("can't read in potential data file %s.\n", filename);
	exit(0);
    }
    printf("opening data files and reading them in.\n");
    numread = fread(MagPot, sizeof(double), NUMELMNTS, readmag);
    if(numread < NUMELMNTS)
	printf("\nWARNING: not all field file read in. \n\n");
    fclose(readmag);

/*  loop over radius  and create data for each one  */

    for(i=1; i<MAXSTEPS+1; i++)
    {

/*  create file to hold computed output from this program  */

	sprintf(filename, "magPot_%02d.dat", i);
	writedat = fopen(filename, "w");
	if (!writedat)
	{
	    printf("can't create electric potential plot file. \n");
	    exit(0);
	}
	idex = i*(i+1)*(i+2)/6;
	across = pi4/i;
	for(k=0; k<=i; k++)
	{
	    kdex = k*(k+1)/2 + idex;
	    phi = across * k;  //  angle from x axis
	    for( j=0; j<=k; j++)
	    {
		theta = across * j;  //  angle from x,y plane
		fprintf(writedat, "%lf	%lf	%lf\n", phi, theta, MagPot[3*(kdex+j)+2]);
	    }
	}
	fclose(writedat);
    }
}