#include "peg_pic.h"
/***************************************************************/  
/* Dumping the current state of the system into a binary file. */
/* Note: the binary file is written in the format used by IBM  */
/* and DEC compatible machines (low bytes followed by high     */
/* bytes). For this Dump() and Restore() call the functions    */
/* sun_read() and sun_write().                                 */

/* In order to have a device independent dmp file you have to
   choose origdmp=origrest=0; Then the dmp file will be in ASCII
   format. In this case size of dmp file can be 2-3 times larger
   than the one written in a binary format.
 */
/* 
 * from 30.7.09 - file also dumps current values of potential into dump file */

char Revision[]={'1','.','1','2'};

/* definition of functions that are compiled only when NOT using xgrafix; otherwise these functions are defined there. */
#ifdef WITH_XGRAFIX

#include "xgrafix/xgrafix.h"

#else

void sun_write(char *, int, int, FILE *);
void sun_read(char *, int, int, FILE *);
void XGWrite(void *, int, int, FILE *, char *);
void XGRead(void *, int, int, FILE *, char *);

#endif


/*
 * Create dump file in binary or ASCII format
 * if origdmp == 3, history of dump files is created; but target filename has to have step inside it (defined in main program)
 */
void Dump( char *filename) {
	int i, j, ftempint;
	float ftemp;
	double lftemp;
	FILE *DMPFile;

	if ((DMPFile = fopen(filename, "w")) == NULL) {
		printf("Dump: open of file %s failed.\n", filename );
		return;
	}

	if( origdmp & 1 ) {
		/* binary dump file */
		XGWrite(Revision, 1, 4, DMPFile, "char");
		ftemp = t;
		XGWrite(&ftemp, 4, 1, DMPFile, "float");
		XGWrite(&oldsigma, 4, 1, DMPFile, "float");
		XGWrite(&extq, 4, 1, DMPFile, "float");
		XGWrite(&extq_1, 4, 1, DMPFile, "float");
		XGWrite(&extq_2, 4, 1, DMPFile, "float");
		XGWrite(&extq_3, 4, 1, DMPFile, "float");
		XGWrite(&nsp, 4, 1, DMPFile, "int");
		XGWrite(jwall, 4, nsp, DMPFile, "float");
		XGWrite(np, 4, nsp, DMPFile, "int");

		for (j=0; j<nsp; j++) {
			for (i=0; i<np[j]; i++) {
				lftemp = x[j][i]/xnc;
				XGWrite(&lftemp, 8, 1, DMPFile, "double");
				ftemp = vscale*vx[j][i];
				XGWrite(&ftemp, 4, 1, DMPFile, "float");
				ftemp = vscale*vy[j][i];
				XGWrite(&ftemp, 4, 1, DMPFile, "float");
				ftemp = vscale*vz[j][i];
				XGWrite(&ftemp, 4, 1, DMPFile, "float");
			}
		}
		
		/* write current values of potential */
		/* first write number of values */
		XGWrite(&nc, 4, 1, DMPFile, "int");
		
		for( i=0 ; i < nc ; i++ ){
			/* write all values */
			XGWrite(&phi[i], 4, 1, DMPFile, "float");
		}
		
	} else {
		/* ASCII dump file */
		lftemp = t;
		fprintf(DMPFile, "%g\n", lftemp);
		ftemp = oldsigma;
		fprintf(DMPFile, "%g\n", ftemp);
		ftemp = extq;
		fprintf(DMPFile, "%g\n", ftemp);
		ftemp =extq_1;
		fprintf(DMPFile, "%g\n", ftemp);
		ftemp =extq_2;
		fprintf(DMPFile, "%g\n", ftemp);
		ftemp =extq_3;
		fprintf(DMPFile, "%g\n", ftemp);
		ftempint= nsp;
		fprintf(DMPFile, "%d\n", ftempint);

		for (j=0; j<nsp; j++) {
			ftemp =jwall[j];
			fprintf(DMPFile, "%g\n", ftemp);
		}

		for (j=0; j<nsp; j++) {
			ftempint=np[j];
			fprintf(DMPFile, "%d\n", ftempint);
		}

		for (j=0; j<nsp; j++) {
			for (i=0; i<np[j]; i++) {
				lftemp = x[j][i];
				fprintf(DMPFile, "%lf\n", lftemp);
				ftemp = vx[j][i];
				fprintf(DMPFile, "%g\n", ftemp);
				ftemp = vy[j][i];
				fprintf(DMPFile, "%g\n", ftemp);
				ftemp = vz[j][i];
				fprintf(DMPFile, "%g\n", ftemp);
			}
		}

		/* write current values of potential */
		/* first write number of values */
		fprintf(DMPFile, "%d\n", nc);
		
		for( i=0 ; i < nc ; i++ ){
			/* write all values */
			fprintf(DMPFile, "%g\n", phi[i]);
		}
	
	}

	fclose(DMPFile);
}
/* END of Dump */


/********************************************

  Function: DumpCoords

  Create some kind of dump file with current values of coordinates (x) and velocities (Vx, Vy, Vz).
  Create separate file for each species

 **********************************************/
void DumpCoords( char *fnbase, char *fnbaseAll, int count) {
	int i, j;
	char tmpfname[255], tmpfname_cum[255];
	FILE *DMPFile, *DMPFile_cum;


	if( origdmp & 2 ) {
		/* write ASCII data */
		for (j=0; j<nsp; j++) {
			/* dump coordinates and velocities for all species */
			if( origdmp > 2 ) {
				/* if more than 2, write separate file for each time step */
				/* create file name */
				sprintf(tmpfname, fnbase, j, count );
				sprintf(tmpfname_cum, fnbaseAll, j );
				
			} else {
				/* write only one file - with each time step, file is overwritten */
				/* create file name */
				sprintf(tmpfname, fnbase, j, 0 );
				sprintf(tmpfname_cum, fnbaseAll, j );
			}

			if ((DMPFile = fopen(tmpfname, "w")) == NULL) {
				printf("DumpCoords: file %s open failed.", tmpfname);
				return;
			}

			if ((DMPFile_cum = fopen(tmpfname_cum, "a")) == NULL) {
				printf("DumpCoords: file %s open failed.", tmpfname_cum);
				return;
			}

			for (i=0; i<np[j]; i++) {
				fprintf( DMPFile, "%d %lf %g %g %g\n", j, (double) length * x[j][i]/nc, vx[j][i]*vscale, vy[j][i]*vscale, vz[j][i]*vscale );
				fprintf( DMPFile_cum, "%d %lf %g %g %g\n", j, (double) length * x[j][i]/nc, vx[j][i]*vscale, vy[j][i]*vscale, vz[j][i]*vscale );
			}

			fclose(DMPFile);
			fclose(DMPFile_cum);
		}
	}

}
/* END of DumpCoords */


/****************************************************************/

void Restore( char *filename) {
	char Rev[5];
	int i, j, ftempint, tmp_nc;
	float ftemp;
	double lftemp;
	FILE *DMPFile;

	if ((DMPFile = fopen(filename, "r+b")) == NULL) {
		puts("Dump: open failed");
		return;
	}
	
	if ( origdmp & 1 ) {
		XGRead(Rev, 1, 4, DMPFile, "char");
		for (i=0; i<4; i++) {
			if (Rev[i]!=Revision[i]) {
				puts("Incompatible dump file version");
				putchar(7);
				exit(1);
			}
		}

		XGRead(&ftemp, 4, 1, DMPFile, "float");
		t = ftemp;
		XGRead(&sigma, 4, 1, DMPFile, "float");
		XGRead(&extq, 4, 1, DMPFile, "float");
		XGRead(&extq_1, 4, 1, DMPFile, "float");
		XGRead(&extq_2, 4, 1, DMPFile, "float");
		XGRead(&extq_3, 4, 1, DMPFile, "float");
		XGRead(&nsp, 4, 1, DMPFile, "int");
		XGRead(jwall, 4, nsp, DMPFile, "float");
		XGRead(np, 4, nsp, DMPFile, "int");

		for (j=0; j<nsp; j++) {
			for (i=0; i<np[j]; i++)	{
				XGRead(&lftemp, 8, 1, DMPFile, "double");
				x[j][i] = lftemp*xnc;
				XGRead(&ftemp, 4, 1, DMPFile, "float");
				vx[j][i]= ftemp/vscale;
				XGRead(&ftemp, 4, 1, DMPFile, "float");
				vy[j][i]= ftemp/vscale;
				XGRead(&ftemp, 4, 1, DMPFile, "float");
				vz[j][i]= ftemp/vscale;
			}
		}

		/* read current values of potential */
		/* first read number of values */
		XGRead(&tmp_nc, 4, 1, DMPFile, "int");
		
		for( i=0 ; i < tmp_nc ; i++ ){
			/* read all values */
			XGRead(&ftemp, 4, 1, DMPFile, "float");
			phi[i] = ftemp;
		}
		
		/* end of if (origdmp) */
	} else {
		fscanf(DMPFile, "%lf\n", &lftemp);
		t=lftemp;
		fscanf(DMPFile, "%g\n", &ftemp);
		sigma=ftemp;
		fscanf(DMPFile, "%g\n", &ftemp);
		extq=ftemp;
		fscanf(DMPFile, "%g\n", &ftemp);
		extq_1=ftemp;
		fscanf(DMPFile, "%g\n", &ftemp);
		extq_2=ftemp;
		fscanf(DMPFile, "%g\n", &ftemp);
		extq_3=ftemp;
		fscanf(DMPFile, "%d\n", &ftempint);
		nsp=ftempint;

		for (j=0; j<nsp; j++) {
			fscanf(DMPFile, "%g\n", &ftemp);
			jwall[j]=ftemp;
		}

		for (j=0; j<nsp; j++) {
			fscanf(DMPFile, "%d\n", &ftempint);
			np[j]=ftempint;
		}

		for (j=0; j<nsp; j++) {
			for (i=0; i<np[j]; i++)	{
				fscanf(DMPFile, "%lf\n", &lftemp);
				x[j][i]=lftemp;
				fscanf(DMPFile, "%g\n", &ftemp);
				vx[j][i]=ftemp;
				fscanf(DMPFile, "%g\n", &ftemp);
				vy[j][i]=ftemp;
				fscanf(DMPFile, "%g\n", &ftemp);
				vz[j][i]=ftemp;
			}
		}

		/* read current values of potential */
		/* first read number of values */
		fscanf(DMPFile, "%d\n", &tmp_nc);
		
		for( i=0 ; i < tmp_nc ; i++ ){
			/* read all values */
			fscanf(DMPFile, "%g\n", &ftemp);
			phi[i] = ftemp;
		}
	
	}   /* end of else */
	fclose(DMPFile);
}
/* END of Restore */

/****** creating file with diagnostics in ASCII format ******/
/*
 * Files created depend on value of origdmp variable from input file (a little bit non logical but what can we do?)
 * for origdmp == 1, only one file is created (no history)
 * for origdmp == 3, one DAT is created and same data written into following files: block1, block2, block3 + all species
 * for origdmp == 5, for whole DAT file, only one file is written (as for origdmp==1), history is made for block1  
 */
void cridatfile(int datfile, char *datfilename, char *datfilenameAll, int currNum) {
	int   i, j, k, isp;
	float vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA;
	FILE  *mydmp, *mydmpBlock;
	char datfname[80], datfname1[80], tmpc[80];


	if( origdmp == 3 ) {
		/* define name for main DAT file (the one that holds all data) */
		sprintf(datfname, datfilenameAll, currNum );
	} else if( origdmp == 1 || origdmp == 5 ){
		/* one single DAT file, no history */
		strcpy(datfname, datfilenameAll);
	}

	/* open file for all the data */
	mydmp = fopen(datfname , "w");

	if (datfile==1) {
		/* open a file for storing only current block of data */
		if( origdmp == 3 || origdmp == 5 ) {
			/* write block 1 data into separate file */
			sprintf(datfname1, datfilename, currNum, "block1" );
			mydmpBlock = fopen(datfname1 , "w");
		} else {
			/* no file opened, define pointer as NULL and use this to check if we can write to file or not */
			mydmpBlock = NULL;
		}

		if (!T_diag) {
			fprintf(mydmp,"t(/dt) = %d, n_grid, V_par (if B>0), V_x, n, T, Qx_tot, Qx_ht, isp=2, ..., Pot., E (if field>0)\n", tstep);
			if( mydmpBlock != NULL) fprintf(mydmpBlock,"t(/dt) = %d, n_grid, V_par (if B>0), V_x, n, T, Qx_tot, Qx_ht, isp=2, ..., Pot., E (if field>0)\n", tstep);
		} else {
			fprintf(mydmp,"t(/dt) = %d, n_grid, V_par (if B>0), V_x, n, T, Qx_tot, Qx_ht, T_par-T_perp, isp=2, ..., Pot., E (if field>0)\n", tstep);
			if( mydmpBlock != NULL) fprintf(mydmpBlock,"t(/dt) = %d, n_grid, V_par (if B>0), V_x, n, T, Qx_tot, Qx_ht, T_par-T_perp, isp=2, ..., Pot., E (if field>0)\n", tstep);
		}

		for(i=0; i<=nc; i++) {                       /* Creating profiles */
			fprintf(mydmp," %d", i);
			if( mydmpBlock != NULL) fprintf(mydmpBlock," %d", i);

			if (bmag>0.0) {
				for (isp=0; isp<nsp; isp++) {
					vv =Vpar_ave_show[isp][i]*dx/dt;
					vv1=Vx_ave_show[isp][i]*dx/dt;
					vv2=srho_ave_show[isp][i]*nc2p/area/dx;
					vv3=Tx_ave_show[isp][i];
					vv4=Qtot_ave_show[isp][i];
					vv5=Qh_ave_show[isp][i];

					if (isp < nsp-1) {
						fprintf(mydmp," %e  %e  %e  %e  %e  %e", vv, vv1, vv2, vv3, vv4, vv5);
						if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e", vv, vv1, vv2, vv3, vv4, vv5);

						if (T_diag) {        
							vv= T_dif_ave_show[isp][i];
							fprintf(mydmp," %e ", vv);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %e ", vv);
						}
					} else {
						if (field) {
							vv6 = phi_ave_show[i];
							vv7= e_ave_show[i];
							fprintf(mydmp," %e  %e  %e  %e  %e  %e", vv, vv1, vv2, vv3, vv4, vv5);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e", vv, vv1, vv2, vv3, vv4, vv5);

							if (T_diag) {
								vv= T_dif_ave_show[isp][i];
								fprintf(mydmp," %e ", vv);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %e ", vv);
							}

							fprintf(mydmp," %e  %e\n", vv6, vv7);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e\n", vv6, vv7);
						} else {
							if (T_diag) {
								vv6= T_dif_ave_show[isp][i];
								fprintf(mydmp," %e  %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6);
							} else {
								fprintf(mydmp," %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5);
							}
						}
					}
				}
			} else {
				for (isp=0; isp<nsp; isp++) {
					vv1=Vx_ave_show[isp][i]*dx/dt;
					vv2=srho_ave_show[isp][i]*nc2p/area/dx;
					vv3=Tx_ave_show[isp][i];
					vv4=Qtot_ave_show[isp][i];
					vv5=Qh_ave_show[isp][i];

					if (isp < nsp-1) {
						fprintf(mydmp,"  %e  %e  %e  %e  %e", vv1, vv2, vv3, vv4, vv5);
						if( mydmpBlock != NULL) fprintf(mydmpBlock,"  %e  %e  %e  %e  %e", vv1, vv2, vv3, vv4, vv5);
					} else {
						if (field) {
							vv6 = phi_ave_show[i];
							vv7 = e_ave_show[i];
							fprintf(mydmp," %e  %e  %e  %e  %e  %e  %e\n", vv1, vv2, vv3, vv4, vv5, vv6, vv7);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e  %e\n", vv1, vv2, vv3, vv4, vv5, vv6, vv7);
						} else {
							fprintf(mydmp," %e  %e  %e  %e  %e\n", vv1, vv2, vv3, vv4, vv5);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e\n", vv1, vv2, vv3, vv4, vv5);
						}
					}
				}
			}
		}  /* end of for(i=0; i<=nc; i++) loop: end of profiles */

		if( mydmpBlock != NULL) {
			fclose(mydmpBlock);
		}

		if( origdmp == 3 ) {
			/* we are writing detailed DAT file + history of it (a lot of data...); open the file */
			/* open a file for storing only current block of data */
			sprintf(datfname1, datfilename, currNum, "block2" );
			mydmpBlock = fopen(datfname1 , "w");
		} else {
			/* no file opened, define pointer as NULL and use this to check if we can write to file or not */
			mydmpBlock = NULL;
		}

		fprintf(mydmp," t/4_pow_x, N(t), LHS_flux, RHS_flux, LHS_Eflux, RHS_Eflux, LHS_T, RHS_T, T(t), T_par_perp(t), Vx(t) (if T_diag>1), gridrefl_left(t), gridrefl_right(t), isp=2, ..., Pot.(1,t), Pot.(2,t), Field_En(t) (if field>0), H_alphaL_1, H_alphaR_1, ..._2, ..._3\n");
		if( mydmpBlock != NULL) fprintf(mydmpBlock," t/4_pow_x, N(t), LHS_flux, RHS_flux, LHS_Eflux, RHS_Eflux, LHS_T, RHS_T, T(t), T_par_perp(t), Vx(t) (if T_diag>1), gridrefl_left(t), gridrefl_right(t), isp=2, ..., Pot.(1,t), Pot.(2,t), Field_En(t) (if field>0), H_alphaL_1, H_alphaR_1, ..._2, ..._3\n");

		//ntm = (hist_hi-1)*interval/nfft1;
		//di = nfft1/interval; 
		//i = 0;
		//for(k=0; k<=ntm; k++, i+=di)        	     /* Time histories */
		for(i=0; i<hist_hi; i++) {
			//fprintf(mydmp," %d", k);
			fprintf(mydmp," %d", i);
			if( mydmpBlock != NULL) fprintf(mydmpBlock," %d", i);

			if (T_diag) {
				for (isp=0; isp<nsp; isp++) {
					j   = np_hist[isp][i];
					vv  = F_hist[isp][i];
					vv1 = F_histRHS[isp][i];
					vv2 = Eflhs_hist[isp][i];
					vv3 = Efrhs_hist[isp][i];
					vv4 = Tlhs_hist[isp][i];
					vv5 = Trhs_hist[isp][i];
					vv6 = Tx_hist[isp][i];
					vv7 = T_dif_hist[isp][i];
					vv8 = Vxt[isp][i]*dx/dt;
					vv9 = gridref_l_hist[isp][i];
					vvA = gridref_r_hist[isp][i];

					if (isp < nsp-1) {
						/* all species but last */
						fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
						if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
					} else {
						/* for last species */
						if (field) {
							fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);

							vv = com_phi_hist[0][i];
							vv1= com_phi_hist[1][i];
							vv2= ese_hist[i];

							if (Nnum) {
								vv3= Halhs_hist[0][i];
								vv4= Harhs_hist[0][i];

								if (Nnum==1) {
									fprintf(mydmp," %f  %f  %e %e %e\n",  vv, vv1, vv2, vv3, vv4);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e %e %e\n",  vv, vv1, vv2, vv3, vv4);
								} else if (Nnum==2) {
									vv5= Halhs_hist[1][i];
									vv6= Harhs_hist[1][i];
									fprintf(mydmp," %f  %f  %e  %e %e %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e  %e %e %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6);
								} else {
									vv5= Halhs_hist[1][i];
									vv6= Harhs_hist[1][i];
									vv7= Halhs_hist[2][i];
									vv8= Harhs_hist[2][i];
									fprintf(mydmp," %f  %f  %e  %e  %e %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e  %e  %e %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8);
								}
							} else {    /* end of if (Nnum)... */
								/* there are no neutrals */
								fprintf(mydmp," %f  %f %e\n",  vv, vv1, vv2);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f %e\n",  vv, vv1, vv2);
							} 
						} else { /* end of if (field) loop */
							if (Nnum) { 
								/* E=0, there are neutrals */
								fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);

								vv2= Halhs_hist[0][i];
								vv3= Harhs_hist[0][i];

								if (Nnum==1) {
									fprintf(mydmp," %e  %e\n",   vv2, vv3);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e\n",   vv2, vv3);
								} else if (Nnum==2) {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									fprintf(mydmp," %e  %e  %e  %e\n",  vv2, vv3, vv4, vv5);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e\n",  vv2, vv3, vv4, vv5);
								} else {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									vv6= Halhs_hist[2][i];
									vv7= Harhs_hist[2][i];
									fprintf(mydmp," %e  %e  %e  %e  %e  %e\n", vv2, vv3, vv4, vv5, vv6, vv7);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e\n", vv2, vv3, vv4, vv5, vv6, vv7);
								}
							} else {    /* end of if (Nnum)... */
								/* E=0, no neutrals */
								fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7, vv8, vv9, vvA);
							}
						}
					}    /* end of else, i.e. after the last particle species */
				}    /* end of for (isp=0; isp<nsp; isp++) loop */
			}  else {   /* end of if (T_diag) loop */
				/* no V(t), T(t), T_par-T_perp diagnostics */
				for (isp=0; isp<nsp; isp++) {
					j   = np_hist[isp][i];
					vv = F_hist[isp][i];
					vv1 = F_histRHS[isp][i];
					vv2 = Eflhs_hist[isp][i];
					vv3 = Efrhs_hist[isp][i];
					vv4 = Tlhs_hist[isp][i];
					vv5 = Trhs_hist[isp][i];
					vv6 = gridref_l_hist[isp][i];
					vv7 = gridref_r_hist[isp][i];

					if (isp < nsp-1) {
						/* all but ast species */
						fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
						if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
					} else {
						/* last particle species */
						if (field) {
							fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
							if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);

							vv = com_phi_hist[0][i];
							vv1= com_phi_hist[1][i];

							if (Nnum) {
								vv2= Halhs_hist[0][i];
								vv3= Harhs_hist[0][i];

								if (Nnum==1) {
									fprintf(mydmp," %f  %f  %e  %e\n",  vv, vv1, vv2, vv3);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e  %e\n",  vv, vv1, vv2, vv3);
								} else if (Nnum==2) {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									fprintf(mydmp," %f  %f  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5);
								} else {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									vv6= Halhs_hist[2][i];
									vv7= Harhs_hist[2][i];
									fprintf(mydmp," %f  %f  %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f  %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
								}
							} else {/* end of if (Nnum)... */
								/* there are no neutrals */
								fprintf(mydmp," %f  %f\n",  vv, vv1);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %f  %f\n",  vv, vv1);
							}
						} else { /* end of if (field) loop */
							if (Nnum) { 
								/* E=0, there are neutrals */
								fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);

								vv2= Halhs_hist[0][i];
								vv3= Harhs_hist[0][i];

								if (Nnum==1) {
									fprintf(mydmp," %e  %e\n",   vv2, vv3);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e\n",   vv2, vv3);
								}  else if (Nnum==2) {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									fprintf(mydmp," %e  %e  %e  %e\n",  vv2, vv3, vv4, vv5);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e\n",  vv2, vv3, vv4, vv5);
								} else {
									vv4= Halhs_hist[1][i];
									vv5= Harhs_hist[1][i];
									vv6= Halhs_hist[2][i];
									vv7= Harhs_hist[2][i];
									fprintf(mydmp," %e  %e  %e  %e  %e  %e\n", vv2, vv3, vv4, vv5, vv6, vv7);
									if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e  %e  %e  %e\n", vv2, vv3, vv4, vv5, vv6, vv7);
								}
							} else { /* end of if (Nnum)... */
								/* E=0, no neutrals */
								fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
								if( mydmpBlock != NULL) fprintf(mydmpBlock," %d  %e  %e  %e  %e  %e  %e  %e  %e ", j, vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
							}
						}	
					}    /* end of else, i.e. after the last particle species */
				}    /* end of for (isp=0; isp<nsp; isp++) loop */
			}    /* end of else */
		}   /* end of for(i=0; i<hist_hi; i++) loop: end of time histories */

		if( mydmpBlock != NULL) {
			fclose(mydmpBlock);
		}

		for (isp=0; isp<nsp; isp++) {                  /* Particle distributions */
			/* open a file for storing only current block of data */
			if( origdmp == 3 ) {
				/* we are writing detailed DAT file + history of it (a lot of data...); open the file */
				sprintf(tmpc, "block3-sp%02da", isp );
				sprintf(datfname1, datfilename, currNum, tmpc );
				mydmpBlock = fopen(datfname1 , "w");
			} else {
				/* no file opened, define pointer as NULL and use this to check if we can write to file or not */
				mydmpBlock = NULL;
			}

			fprintf(mydmp," isp= %d\n", isp);
			if( mydmpBlock != NULL) fprintf(mydmpBlock," isp= %d\n", isp);
			fprintf(mydmp," E, f(E) LHS Wall, f(E) RHS Wall, f1(E), f2(E), ...\n");
			if( mydmpBlock != NULL) fprintf(mydmpBlock," E, f(E) LHS Wall, f(E) RHS Wall, f1(E), f2(E), ...\n");

			for(i=0; i<nbin[isp]; i++) {
				vv = e_array[isp][i]*Escale[isp];
				vv1 = fe[isp][i];
				vv2 = feR[isp][i];

				if (!NDP) {
					fprintf(mydmp," %e  %e  %e\n", vv, vv1, vv2);
					if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e\n", vv, vv1, vv2);
				} else {
					fprintf(mydmp," %e  %e  %e", vv, vv1, vv2);
					if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e", vv, vv1, vv2);

					for (k=0; k<NDP-1; k++) {
						vv = fe_mid_show[isp][k][i];
						fprintf(mydmp," %e", vv);
						if( mydmpBlock != NULL) fprintf(mydmpBlock," %e", vv);
					}

					vv = fe_mid_show[isp][NDP-1][i];
					fprintf(mydmp," %e\n", vv);
					if( mydmpBlock != NULL) fprintf(mydmpBlock," %e\n", vv);
				}
			}

			if( mydmpBlock != NULL) {
				fclose(mydmpBlock);
			}

			if (NDP) {
				/* open a file for storing only current block of data */
				if( origdmp == 3 ) {
					/* we are writing detailed DAT file + history of it (a lot of data...); open the file */
					sprintf(tmpc, "block3-sp%02db", isp );
					sprintf(datfname1, datfilename, currNum, tmpc );
					mydmpBlock = fopen(datfname1 , "w");
				} else {
					/* no file opened, define pointer as NULL and use this to check if we can write to file or not */
					mydmpBlock = NULL;
				}

				fprintf(mydmp," v, f1(v), f2(v), ...\n");
				if( mydmpBlock != NULL) fprintf(mydmpBlock," v, f1(v), f2(v), ...\n");

				for(i=0; i<nbin_fvx[isp]; i++) {
					vv = fvx_array[isp][i]*vscale;
					fprintf(mydmp," %e", vv);
					if( mydmpBlock != NULL) fprintf(mydmpBlock," %e", vv);

					for (k=0; k<NDP-1; k++) {
						vv = f_vpar_show[isp][k][i];
						fprintf(mydmp," %e", vv);
						if( mydmpBlock != NULL) fprintf(mydmpBlock," %e", vv);
					}

					vv = f_vpar_show[isp][NDP-1][i];
					fprintf(mydmp," %e\n", vv);
					if( mydmpBlock != NULL) fprintf(mydmpBlock," %e\n", vv);
				}

				if( mydmpBlock != NULL) {
					fclose(mydmpBlock);
				}
			}

		}/* end of particle distribution part */

		if (field) {                       /* Oscillation frequency diagnostics */
			/* open a file for storing only current block of data */
			if( origdmp == 3 ) {
				/* we are writing detailed DAT file + history of it (a lot of data...); open the file */
				sprintf(datfname1, datfilename, currNum, "block4" );
				mydmpBlock = fopen(datfname1 , "w");
			} else {
				/* no file opened, define pointer as NULL and use this to check if we can write to file or not */
				mydmpBlock = NULL;
			}

			fprintf(mydmp," Freq., Amp_LHS_Pot, Amp_RHS_Pot\n");
			if( mydmpBlock != NULL) fprintf(mydmpBlock," Freq., Amp_LHS_Pot, Amp_RHS_Pot\n");

			for(i=0; i<freq_hi; i++) {
				vv  = f_array[i];
				vv1 = phi_fft[i];
				vv2= mphi_fft[i];
				fprintf(mydmp," %e  %e  %e\n", vv, vv1, vv2);
				if( mydmpBlock != NULL) fprintf(mydmpBlock," %e  %e  %e\n", vv, vv1, vv2);
			}

			if( mydmpBlock != NULL) {
				fclose(mydmpBlock);
			}
		}
		/* end of datfile=1 */

	} else if (datfile==2) {
		fprintf(mydmp,"t(/dt) = %d\n", tstep);
		fprintf(mydmp," t/4_pow_x, N(t), LHS_flux, RHS_flux, LHS_Eflux, RHS_Eflux, LHS_T, RHS_T, isp=2, ..., Pot.(1,t), Pot.(2,t) (if field>0), H_alphaL_1, H_alphaR_1, ..._2, ..._3\n");

		//ntm = (hist_hi-1)*interval/nfft1;
		//di = nfft1/interval;
		//i = 0;
		//for(k=0; k<=ntm; k++, i+=di)                   /* Time histories */
		for(i=0; i<hist_hi; i++) {
			// fprintf(mydmp," %d", k);
			fprintf(mydmp," %d", i);

			for (isp=0; isp<nsp; isp++) {
				j   = np_hist[isp][i];
				vv = F_hist[isp][i];
				vv1 = F_histRHS[isp][i];
				vv2 = Eflhs_hist[isp][i];
				vv3 = Efrhs_hist[isp][i];
				vv4 = Tlhs_hist[isp][i];
				vv5 = Trhs_hist[isp][i];

				if (isp < nsp-1) {
					/* for all species, but last */
					fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5);

				} else {
					/* for last one */
					if (field) {
						fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5);

						vv = com_phi_hist[0][i];
						vv1= com_phi_hist[1][i];

						if (Nnum) {
							vv2= Halhs_hist[0][i];
							vv3= Harhs_hist[0][i];

							if (Nnum==1) {
								fprintf(mydmp," %f  %f  %e  %e\n",  vv, vv1, vv2, vv3);
							} else if (Nnum==2) {
								vv4= Halhs_hist[1][i];
								vv5= Harhs_hist[1][i];
								fprintf(mydmp," %f  %f  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5);
							} else {
								vv4= Halhs_hist[1][i];
								vv5= Harhs_hist[1][i];
								vv6= Halhs_hist[2][i];
								vv7= Harhs_hist[2][i];
								fprintf(mydmp," %f  %f  %e  %e  %e  %e  %e  %e\n", vv, vv1, vv2, vv3, vv4, vv5, vv6, vv7);
							}
						} else {    /* end of if (Nnum)... */
							/* there are no neutrals */
							fprintf(mydmp," %f  %f\n",  vv, vv1);
						}
						/* end of if (field) loop */
					} else {
						if (Nnum) {
							/* E=0, there are neutrals */
							fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e", j, vv, vv1, vv2, vv3, vv4, vv5);

							vv2= Halhs_hist[0][i];
							vv3= Harhs_hist[0][i];

							if (Nnum==1) {
								fprintf(mydmp," %e  %e\n",   vv2, vv3);
							} else if (Nnum==2) {
								vv4= Halhs_hist[1][i];
								vv5= Harhs_hist[1][i];
								fprintf(mydmp," %e  %e  %e  %e\n",  vv2, vv3, vv4, vv5);
							} else {
								vv4= Halhs_hist[1][i];
								vv5= Harhs_hist[1][i];
								vv6= Halhs_hist[2][i];
								vv7= Harhs_hist[2][i];
								fprintf(mydmp," %e  %e  %e  %e  %e  %e\n", vv2, vv3, vv4, vv5, vv6, vv7);
							}
						} else {    /* end of if (Nnum)... */
							/* E=0, no neutrals */
							fprintf(mydmp," %d  %e  %e  %e  %e  %e  %e\n", j, vv, vv1, vv2, vv3, vv4, vv5);
						}    
					}
				}    /* end of else, i.e. after the last particle species */
			}    /* end of for (isp=0; isp<nsp; isp++) loop */
		}    /* end of time histories */
	} else {    /* end of datfile = 2 */
		/* Datfile 3 */
		fprintf(mydmp,"t(/dt) = %d\n", tstep);
		fprintf(mydmp," t/4_pow_x, N(t), LHS_flux, RHS_flux, isp=2, ..., Pot.(1,t), Pot.(2,t) (if field>0)\n");

		//ntm = (hist_hi-1)*interval/nfft1;
		//di = nfft1/interval;
		//i = 0;
		//for(k=0; k<=ntm; k++, i+=di)                   /* Time histories */
		for(i=0; i<hist_hi; i++) {
			//fprintf(mydmp," %d", k);
			fprintf(mydmp," %d", i);

			for (isp=0; isp<nsp; isp++) {
				j   = np_hist[isp][i];
				vv = F_hist[isp][i];
				vv1 = F_histRHS[isp][i];

				if (isp < nsp-1) {
					fprintf(mydmp," %d  %e  %e", j, vv, vv1);

				} else {
					if (field) {
						fprintf(mydmp," %d  %e  %e", j, vv, vv1);

						vv = com_phi_hist[0][i];
						vv1= com_phi_hist[1][i];
						fprintf(mydmp," %f  %f\n",  vv, vv1);
					} else {
						fprintf(mydmp," %d  %e  %e\n", j, vv, vv1);
					}
				}
			}
		}
	} /* end of else (datfile 3) */

	fclose(mydmp);
}       
/* END of cridatfile */

/****************************************************/

/*
 * functions from xgrafix/xgrafix.c
 * 
 * They are compiled only, when XGrafix is NOT compiled in - just to be able to preserve same function from writing/reading dump and dat files
 * 
 */

#ifndef WITH_XGRAFIX

void XGRead(void *ptr, int size, int nitems, FILE *stream, char *type) {
#ifdef IBM
  sun_read(ptr, size, nitems, stream);
#endif

#ifdef DEC
  fread(ptr, size, nitems, stream);
#endif

#ifdef SUN
  sun_read(ptr, size, nitems, stream);
#endif

#ifdef HP
  sun_read(ptr, size, nitems, stream);
#endif

#ifdef UNICOS
  if(!strcmp(type, "float"))     unicos_read_float(ptr, nitems, stream);
  else if(!strcmp(type, "int"))  unicos_read_int(ptr, nitems, stream);
  else if(!strcmp(type, "char")) fread(ptr, size, nitems, stream);
  else                           puts("XGRead: unrecognized type");
#endif
}

/******************************************************/

void XGWrite(void *ptr, int size, int nitems, FILE *stream, char *type) /* ddd*/
{
#ifdef IBM
  sun_write(ptr, size, nitems, stream);
#endif
  
#ifdef DEC
  fwrite(ptr, size, nitems, stream);
#endif
  
#ifdef SUN
  sun_write(ptr, size, nitems, stream);
#endif
  
#ifdef HP
  sun_write(ptr, size, nitems, stream);
#endif

#ifdef UNICOS
  if(!strcmp(type, "float"))     unicos_write_float(ptr, nitems, stream);
  else if(!strcmp(type, "int"))  unicos_write_int(ptr, nitems, stream);
  else if(!strcmp(type, "char")) fwrite(ptr, size, nitems, stream); 
  else                           puts("XGWrite: unrecognized type");
#endif
}

/******************************************************/

void sun_read(char *ptr, int size, int nitems, FILE *stream) /* ddd */
{
  int i, j;
  
  for(j= 0; j <nitems; j++)
    for(i= size-1; i >=0; i--)
      fread(ptr+j*size+i, 1, 1, stream);
}

/******************************************************/

void sun_write(char *ptr, int size, int nitems, FILE *stream) /* ddd */
{
  int i, j;
  
  for(j= 0; j <nitems; j++)
    for(i= size-1; i >=0; i--)
      fwrite(ptr+j*size+i, 1, 1, stream);
}

/******************************************************/

#ifdef UNICOS
void unicos_read_float(float *ptr, int nitems, FILE *stream) /* ddd */
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;

  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    fread(&d,1,1,stream);
    fread(&c,1,1,stream);
    fread(&b,1,1,stream);
    fread(&a,1,1,stream);

    if( a == 0 && b == 0 && c == 0 && d == 0) *iptr = 0;
    else 
    {
      t = ((a << 1) | (b >> 7));
      *iptr = ((a  >> 7) << 63) | ((t -126 + 040000) << 48);
      *iptr |= (b | 0x80) << 40;  
      *iptr |= c << 32;
      *iptr |= d << 24;
    }
  }
}
#endif

/******************************************************/

#ifdef UNICOS
void unicos_read_int(int *ptr, int nitems, FILE *stream)   /* ddd */
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;

  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    fread(&a,1,1,stream);
    fread(&b,1,1,stream);
    fread(&c,1,1,stream);
    fread(&d,1,1,stream);

    *iptr = d << 24 | c << 16 | b << 8 | a;
    *iptr |= (d >> 7) ? 0xffffffff00000000 : 0x0;
  }
}
#endif

/******************************************************/

#ifdef UNICOS
void unicos_write_float(float *ptr, int nitems, FILE *stream)   /* ddd */
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;
  
  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);

    t = ((*iptr >> 48) & 0xff) - 040000 + 126;
    
    a = ((*iptr & (1 << 63)) ? 1 << 7 : 0) | t >> 1; 
    b = (t & 0x1) << 7 | ((*iptr >> 40) & 0x7f);
    c = *iptr >> 32;
    d = *iptr >> 24;
    
    fwrite(&d,1,1,stream);
    fwrite(&c,1,1,stream);
    fwrite(&b,1,1,stream);
    fwrite(&a,1,1,stream);
  }
}
#endif

/******************************************************/

#ifdef UNICOS
void unicos_write_int(int *ptr, int nitems, FILE *stream)   /* ddd */
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;
  
  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    a = *iptr & 0xff;
    b = *iptr >> 8 & 0xff;
    c = *iptr >> 16 & 0xff;
    d = *iptr >> 24 & 0xff;
    fwrite(&a,1,1,stream);
    fwrite(&b,1,1,stream);
    fwrite(&c,1,1,stream);
    fwrite(&d,1,1,stream);
  }
}
#endif

#endif
