#include "peg_pic.h"

/* Creating the file for "movie". */
/* "tchar" is the name of the criated file: <tchar>.movie.txt; */
/* "k" is the time interval (in time steps) between snap shots */
/* "dk" is the time interval for time averaging of profiles.   */ 
/* The routine calculates time when the averaging has to be    */
/* started and finished, calculates dk time step averaged n, T */
/* and Q profiles and writes them into the <tchar>.movie.txt   */
/* file.                                                       */ 


void movf(char *tchar, int k, int dk) {
	int           i, j, s;
	float         vv, vv1, vv2, vv3, vv4, vv5;
	static int    initflag=1, mt=0, nstp=0;
	static float  ***outt;
	static char   nm[80];
	FILE          *mydmp;

	nstp +=1;

	if ((!qtp) && (nstp==k-1)) {
		qtp=3;         /* Flag for calculation of n, T and Q profiles */
	}

	if (initflag) {                    /* Initialization */
		strcpy(nm, tchar);
		strcat(nm, ".movie.txt");

		/* Create array for calculation of averaged profiles */
		outt= (float ***)malloc(6*sizeof(float **));

		for(i=0; i<6; i++) {
			outt[i]= (float **)malloc(nsp*sizeof(float *));

			for(j=0; j<nsp; j++) {
				outt[i][j]= (float *)malloc(ng*sizeof(float));
			}
		}

		for(i=0; i<6; i++) {
			for(j=0; j<nsp; j++) {
				for(s=0; s<ng; s++) {
					outt[i][j][s]= 0;
				}
			}
		}

		// why    mydmp = fopen(nm , "w");     /* First delete the existing file */
		//         fclose(mydmp);

		nstp= k;                /* For calcilation of "initial" profiles */
		initflag=0;
	}


	if (nstp==k) {
		mt=dk;       /* Start counting for averaging of profiles */
	}
	
	if (mt) {
		for(i=0; i<nsp; i++) {
			for(j=0; j<ng; j++) {
				outt[0][i][j] += srho[i][j];
				outt[1][i][j] += Vpar[i][j];
				outt[2][i][j] += Tx[i][j];
				outt[3][i][j] += Qtot[i][j]; 	
				outt[4][i][j] += Qh[i][j];
			}
		}

		for(j=0; j<ng; j++) {
			// potencial is the same for all particles so just store value into first cell
			outt[5][0][j] += phi_ave_show[j];
		}

		if (mt==1) {                     /* The last step for averaging */
			mydmp = fopen(nm , "a");

			for(i=0; i<ng; i++) {
				fprintf(mydmp,"%d", i);

				for(j=0; j<nsp; j++) {
					vv = outt[0][j][i]*nc2p/area/dx/dk;
					vv1= outt[1][j][i]/dk;
					vv2= outt[2][j][i]/dk;
					vv3= outt[3][j][i]/dk;
					vv4= outt[4][j][i]/dk;

					fprintf(mydmp,"  %e  %e  %e  %e  %e", vv, vv1, vv2, vv3, vv4);
				}

				vv5= outt[5][0][i]/dk;
				fprintf(mydmp,"  %e  %e\n", vv5, phi_ave_show[i]);
			}

			fclose(mydmp);

			for(i=0; i<5; i++) {
				for(j=0; j<nsp; j++) {
					for(s=0; s<ng; s++) {
						outt[i][j][s]= 0;
					}
				}
			}

			if (qtp==3) {
				qtp=0;         /* No profiles until next "movie" step */
			}
			
			nstp= 0;                  /* Initialize the counter */
		}

		mt--;
	} /* end of if (mt) */

} /* END of movef */

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