/************
* sphere.c 
*          Steven L. Jacques, Oct. 1, 2002
* Plots p(z,x) at given time t
* for a series of time points.
* This program is written to simulate uniform
* energy deposition in a sphere at time zero.
* It can be modified by adjusting the specification of W.
* Yields files
*	pzx##.dat for ##=0:40 timepoints, zzx.dat, and Wzx.dat
* which are graphed using another application,
* for example, MATLAB.
* Each timepoint takes about 8 min to calculate.
*****/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "array.h"
#include "nrutil.h"
#include <time.h>

#define	Nz      	200
#define	Nx      	200

/*************************************************************/	
/*************** MAIN ****************************************/	
int main(void) {
	/* constants */
	double pi   = 3.14159265358979;
	double cs   = 1500;           /* [m/s] */
	double rho  = 1000;           /* [kg] */
	double Cp   = 4180;           /* [J/(kg deg)] */
	double G    = 0.12;           /* [-] */
	/* double mua  = 0.5e4;            [m^-1] , 1e4 m^-1 = 100 cm^-1*/
	double Wo   = 1e6;            /* [J/m^3 */
	
	/* variables */
	double dz;           /* [m], bin size */
	double r, r2, x, z, zc, dr, dx, dt, beta, t, radius;
	double ur, mvp1, mvp2, summ1, summ2;
	double start_time, finish_time; /* for clock() */
	double dtime = 0;
	double dtt;
	double W, rad2, xp, yp, zp, th, ph, dth, dph;
	time_t now;
	char   name[30];
	
	/* dummy variables */
	double u;
	int    i, j, k, m, it;
	int    Nth = 100;
	int    Nph = 100;
	FILE*	target;     /* point to output file */

	/* arrays */
	double **p;
	p = dmatrix(1, Nz, 1, Nx);
	

printf("begin\n");
/*********** BEGIN ************/

now = time(NULL);
printf("%s\n", ctime(&now));

zc     = 0.0;            /* [m] center of sphere of deposition */
radius = 2e-3;           /* radius of deposition sphere */
dz     = 4*radius/Nz;
dx     = dz;             /* [m] */
dt     = dz/cs;          /* [s] */
dr     = cs*dt;          /* [m] */
beta   = G*Cp/cs/cs;     /* [deg^-1] */
rad2   = radius*radius;  
dtt    = dz*Nz/cs/40;    /* center reaches edge of map */
dth    = pi/Nth;
dph    = 2*pi/Nph;

for (it=0; it<=40; it++) {
	printf("time # %d\n",it);
	t      = dtt*it; /* <------- SPECIFY TIME of snapshot */
	r      = cs*t;
	r2     = cs*(t + dt);
	
/******
* Cycle thru p(z,x) = p[i][j]
*******/
start_time = clock();
for (j=1; j<=Nx; j++) {/* Cycle through x */
	x = dx*(j-0.5);

	if (j>1) {
		finish_time = clock();
		if (j <= 5) dtime = (double)(finish_time - start_time)/CLOCKS_PER_SEC/60*Nx/(j-1)/Nx;
		/*printf("%d\tWill finish in = %5.2f min\n", j, dtime*(Nx-j));
		*/
		}
		
	for (i=1; i<=Nz; i++) {/* Cycle through z */
		z = dz*(i-0.5);
		
		mvp1 = 0.0;
		mvp2 = 0.0;
		for (k=1; k<=Nth; k++) {
			th = (k-0.5)*dth;
			
			summ1 = 0.0;
			summ2 = 0.0;
			for (m=1; m<=Nph; m++) {
				ph = (m - 0.5)*dph;
				
				/* first time point */
				ur = r*sin(th);
				xp = ur*cos(ph) + x;
				yp = ur*sin(ph);
				zp = z - r*cos(th);
				W  = (xp*xp + yp*yp + (zp - zc)*(zp - zc) <= rad2);
				summ1 += W;

				/* second time point */
				ur = r2*sin(th);
				xp = ur*cos(ph) + x;
				yp = ur*sin(ph);
				zp = z - r2*cos(th); /* uz source-to-surface distance */
				W  = (xp*xp + yp*yp + (zp - zc)*(zp - zc) <= rad2);
				summ2 += W;
				
				} /* end m, phi cycle */

			u = sin(th);
			mvp1 += summ1*u;
			mvp2 += summ2*u;
			} /* end k, theta cycle */
			
		mvp1 *= Wo*r*dph*dth*dr*beta/(4*pi*rho*Cp);
		mvp2 *= Wo*r2*dph*dth*dr*beta/(4*pi*rho*Cp);
		p[i][j] = rho*(mvp2 - mvp1)/dt;
		
		} /* end i, z cycle*/

	} /* end j, x cycle */

if (it==1){
	target = fopen("zzx.dat", "w");
	for (i=1; i<=Nz; i++)
		fprintf(target, "%5.4e\n", (i-0.5)*dz);
	fclose(target);
}

sprintf(name, "pzx%d.dat", it);
target = fopen(name, "w");
for (i=1; i<=Nz; i++) {
	for (j=1; j<=Nx; j++) {
		fprintf(target, "%5.4e", p[i][j]);
		if (j<Nx) fprintf(target, "\t");
		}
	fprintf(target, "\n");
	}
fclose(target);

} /***** end it *****/


/*********** END **************/

free_dmatrix(p, 1, Nz, 1, Nx);

finish_time = clock();
printf("--------------------------------------------------------------------\n");
printf("Elapsed Time = %5.2f min\n", (double)(finish_time - start_time)/CLOCKS_PER_SEC/60);
now = time(NULL);
printf("%s\n", ctime(&now));

return(0);
}

/*************** end MAIN ************************************/	
/*************************************************************/