/*##########################################################
# description:	put some data in a 2d field and print it
#		in a two different 3D formats
#########################################################*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>

#define PI 3.1415926535897932

void **alloc2d(int, int, int);
void free2d(void **, int, int);

int main() {

	float	**XZfield;
	int	i, j;
	int	Nx, Nz;
	float	Xmax=10., Zmax=10.;
	float	dX=0.1, dZ=0.1;
	FILE	*fd1, *fd2; 

	Nx = (int) Xmax/dX;
	Nz = (int) Zmax/dZ;

	// allocate 2D matrix for the data
	XZfield	= (float **) alloc2d( Nx, Nz, sizeof(float) );

	// make sure allocation was ok
	if ( XZfield == NULL ) {
		printf("error allocating memory\n");
		exit(1);
	}

	// put zeros all around the edges of the grid
	for (i=0; i<Nx; ++i ) {
		XZfield[i][0] = 0.;
		XZfield[i][Nz-1] = 0.;
	}
	for (j=0; j<Nz; ++j ) {
		XZfield[0][j] = 0.;
		XZfield[Nx-1][j] = 0.;
	}

	// put some data in the grid
	for (i=1; i<Nx-1; ++i ) {
		for (j=1; j<Nz-1; ++j) {
			XZfield[i][j] = sin(0.25*PI*j*dZ) + sin(0.25*PI*i*dX) ;
		}
	}

	// print in x y z style
        if ( (fd1 = fopen("myfield.dgrid3d", "w")) == NULL ) {
                printf("error opening %s\n", "myfield.dgrid3d" );
                exit(1);
        }
	for (i=0; i<Nx; i++)
                for (j=0; j<Nz; j++)
                        fprintf(fd1, "%f %f %f\n", i*dX, j*dZ, XZfield[i][j]);
        fclose(fd1);

	// print in matrix style
        if ( (fd2 = fopen("myfield.pm3d", "w")) == NULL ) {
                printf("error opening %s\n", "myfield.pm3d" );
                exit(1);
        }
	for (i=0; i<Nx; i++) {
                for (j=0; j<Nz; j++)
                        fprintf(fd2, "%f ", XZfield[i][j]);
		fprintf(fd2, "\n");
	}
        fclose(fd2);

	free2d((void **)XZfield, Nx, Nz);
}


void **alloc2d( int d1, int d2, int size) {
	int i,j;
	void **p;
	p = (void**)malloc(d1*sizeof(void*));

	if ( p == NULL ) return NULL;
	for ( i=0; i<d1; ++i ) {
		p[i] = malloc(d2*size);
		if ( p[i] == NULL ) return NULL;
	}
	return p;
}
void free2d( void **p, int d1, int d2 ) {
	int i,j;

	if ( p != NULL ) {
		for ( i=0; i<d1; ++i ) {
			if ( p[i] != NULL )
				free(p[i]);
		}
		free(p);
	}
}